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ABSTRACT 


Shaffer,  M.J.  and  W.E.  Larson,  eds.   1987 
NTRM,  a  Soil-Crop  Simulation  Model  for 
Nitrogen,  Tillage,  and  Crop-Residue 
Management.   U.S.  Department  of 
Agriculture  Conservation  Research  Report 
34-1,  103  pp. 


A  large,  broad-based  computer  simulation 
model  was  developed  for  tillage,  crop- 
residue,  and  nitrogen  fertilizer 
management.   This  NTRM  (Nitrogen-Tillage  - 
Residue  Management)  model  was  designed 
mainly  to  provide  management  assistance 
at  the  farm,  extension,  and  engineering 
levels.   Other  basic  applications  in 
these  areas,  as  well  as  applications  in 
research  and  teaching,  are  also  possible. 
The  NTRM  model  simulates  physical, 
chemical,  and  biological  processes  in  the 
soil-water-crop  continuum  using 
integrated  submodels  for  soil 
temperature,  soil  carbon  and  nitrogen 


transformations,  unsaturated  flow  of 
water,  crop  and  root  growth,  evaporation 
and  transpiration,  tillage,  interception 
and  infiltration,  chemical  equilibria 
processes,  solute  transport,  and  crop 
residues.   The  overall  structure  of  the 
model  is  presented  along  with  submodel 
interactions  and  detailed  technical 
descriptions  of  each  submodel. 
Validation  and  verification  information 
are  provided  for  each  submodel  and  the 
overall  model.   Chapters  are  included 
that  discuss  management  applications  of 
the  model  and  simplified 
versions/configurations  that  may  be 
appropriate  for  specialized  model  uses. 
Future  research  needs  in  this  field  are 
outlined. 

Key  words :   agricultural  management,  crop 
residues,  crop  yield,  infiltration, 
models,  nitrogen,  roots,  soil,  solutes, 
temperature,  tillage,  water 
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Chapter  1 
INTRODUCTION 

M.J.  Shaffer  and  W.E.  Larson1 


In  recent  years ,  a  need  has  been 
identified  to  develop  integrated 
guidelines  for  the  management  of  tillage 
practices,  crop  residues,  and  nitrogen 
fertilizers  in  U.S.  agriculture.   In 
1979,  the  U.S.  Department  of  Agriculture, 
Agricultural  Research  Service  (ARS)  began 
a  national  research  program  to  develop 
these  guidelines.   Several  ARS  research 
units  and  cooperating  universities  have 
been  involved  in  this  work.   The 
development  of  a  comprehensive  computer 
simulation  model  system  for 
nitrogen- tillage -residue  management 
applications  at  various  levels  was  a 
prime  objective  of  the  overall  effort. 
Known  as  the  Nitrogen-Tillage-Residue 
Management  (NTRM)  model.it  was  developed 
by  the  ARS  Soil  and  Water  Management 
Research  Unit  in  St.  Paul,  MN,  with 
assistance  from  the  University  of 
Minnesota,  the  ARS  research  unit  at 
Morris,  MN,  and  other  cooperating 
universities  and  ARS  research  units.   The 
NTRM  model  is  designed  to  have  management 
and  other  basic  applications  at  the 
farmer,  extension,  research,  and 
engineering  levels.   In  addition, 
teaching  of  soil -crop  management  and 
related  applications  are  highly  possible. 


range  of  soil-water-crop  models  can  be 
found  in  the  literature.   The  amount  of 
detail  likely  to  be  found  in  various 
simulation  models  (independent  of 
scientific  field)  is  illustrated  in 
figure  1.1.   Theorists  tend  to  construct 
the  more  detailed  models,  while 
practitioners  favor  the  single-equation 
or  back-of- the -envelope  approaches. 
Models  that  contain  the  most  detail  tend 
to  be  very  flexible  in  the  number  and 
types  of  applications  that  can  be  made, 
provided  anyone  other  than  the  author  can 
understand  them.   On  the  other  hand,  the 
simplest  models  are  easy  to  apply  and 
understand  but  are  limited  in  scope  and 
have  a  narrow  range  of  applicability 
(although  their  wide  usage  often 
disguises  these  facts). 

The  most  successful  models,  however 
(defined  as  those  that  are  used 
frequently  yet  have  good  credibility  and 
applicability) ,  lie  somewhere  in  the 
middle  near  the  maximum  point  on  the 
solid  curve  in  figure  1.1. 

The  NTRM  model  with  all  its  options  in 
use  probably  falls  slightly  toward  the 
theory  side  of  the  optimum.   Less 


Computer  simulation  of  the 
soil-water-crop  continuum  is  a 
multidisciplinary  science  with  its 
beginnings  dating  back  to  the  early 
1960's.   The  degree  of  sophistication  and 
detail  contained  in  these  models  has 
varied  with  the  particular  author,  the 
application,  and  the  development  of 
computer  hardware  and  software.   A  wide 

1  Shaffer  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108,  and 
Larson  is  the  head,  Soil 
Science  Department,  University 
of  Minnesota,  St.  Paul,  MN 
55108. 
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Figure  1.1. 

Probability  of  having  the  best  model  for 

a  given  application. 


vigorous  configurations  of  the  model 
system  can  be  used  by  practitioners.   It 
should  be  noted,  however,  that  this  model 
is  not  an  "all-out"  one,  and  considerable 
effort  has  already  been  made  to  select 
middle-of-the-road  approaches  for  the 
submodels . 

The  NTRM  model  (fig.  1.2)  was  developed 
over  a  period  of  about  4  years  by  using 
and  adapting  some  existing  submodels  and 
by  developing  certain  others.   Working  as 
a  team,  soil  physicists,  soil  chemists, 
hydrologists ,  plant  physiologists,  soil 
microbiologists,  engineers,  and  computer 
scientists  contributed  to  the  development 
of  the  overall  model. 

The  primary  objective  was  to  develop  a 
comprehensive,  physically  significant 
model  of  the  soil  environment  and  its 
effects  on  crop  growth.   Emphasis  was 
placed  on  the  tillage,  nitrogen,  and 
residue  aspects  of  the  system.   The  scope 
of  the  model  was  kept  wide  enough  to 
allow  applications  in  the  areas  of 
research,  crop  and  soil  management,  and 
teaching. 


CHEMICAL 
EQUILIBRIA 
SUBMODEL 


INFILTRATION/ 

SURFACE 

RUNOFF 

SUBMODEL 


CALLS  OTHER 
ROUTINES  ON  DAILY 
AND  0.1  DAY  BASIS 


CALLS  MOST  MAJOR 

COMPUTATIONAL 

ROUTINES 


ROOT 

GROWTH 
SUBMODEL 


3 


CROP  NUTRIENT 

UPTAKE 

SUBMODEL 


NITROGEN 

TRANSFORMATION 
SUBMODELS! 


SALT 
TRANSPORT 
SUBMODEL 


*  t  t  ' ' 


OUTPUT 

TO  DISC  FOR 

PLOTTING  &  LIST 


Figure  1.2. 

Nitrogen- Tillage -Residue -Management 

(NTRM)  Model. 


A  cross  section  of  the  crop-soil-water 
continuum  and  aquifer  being  simulated 
with  the  NTRM  model  is  shown  in  figure 
1.3.   Note  the  basic  physical,  chemical, 
and  biological  processes  that  are  being 
modeled.   The  aquifer  portions  of  the 
model  system  are  a  separate  set  of 
submodels  that  require  input  from  the 
unsaturated- zone  simulators.   Likewise, 
the  soil  temperature  submodel  is  a  stand- 
alone program  that  generates  inputs  for 
the  overall  model. 

The  "Detailed  Return  Flow  Salinity  and 
Nutrient  Simulation  Model"  (fig  1.4) 
reported  by  Shaffer  et  al .  (1977),  was 
used  as  a  core  for  the  NTRM  model.   This 
irrigation  return  flow  model  already 
contained  relatively  comprehensive 
submodels  for  waterflow,  nutrient 
transformations,  salt  reactions,  and  salt 
transport.   One  current  version  of  this 
model,  described  in  Shaffer  and  Gupta 
(1981) ,  has  been  used  to  simulate  solute 
reactions  and  transport  in  land  treatment 
of  municipal  wastewater  effluents. 

The  original  core  model  has  a  history 
dating  back  to  the  early  and  middle 
1960 's  when  G.R.  Dutt  of  the  University 
of  Arizona  and  K.  Tanj i  of  the  University 
of  California,  Davis,  began  developing 
the  first  submodels,  primarily  in  the 
area  of  soil  chemistry  associated  with 
irrigation  in  the  Western  United  States 
(Dutt  1961,  Dutt  and  Tanj i  1962).   By  the 
early  1970' s,  the  modeling  work  had  been 
extended  to  include  unsaturated  flow  and 
nitrogen  transformations  under  a  growing 
crop  (Shaffer  et  al.  1969,  Dutt  et  al . 
1972,  Shaffer  and  Dutt  1973).   The  model 
was  also  extended  by  personnel  at  the 
U.S.  Bureau  of  Reclamation  to  include 
waterflow  and  chemical  processes  in  the 
aquifer  below  irrigated  areas  (Ribbens 
1970,  Shaffer  and  Ribbens  1974,  Shaffer 
et  al.  1977).   In  the  late  1970's,  the 
model  was  adapted  by  the  ARS  group  in  St. 
Paul  to  simulate  nutrient 
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transformations,  solute  chemistry,  and 
waterflow  processes  during  land  treatment 
of  sewage  effluents  (Gupta  et  al .  1978). 

In  late  1979,  work  began  on  the  initial 
development  of  the  NTRM  model  (Shaffer 
1980) .   The  Nebraska  corn  growth  model 
CORNGRO  (Childs  et  al.  1977)  was 
incorporated  with  modifications  as  a 
submodel  for  crop  growth.   Submodels  for 
additional  crops  will  be  added  as  the 
NTRM  model  system  is  further  developed. 
A  dynamic  root  growth  model  (Shaffer  and 
Clapp  1982)  (see  also  Chapter  10  of  this 
handbook)  was  developed  from  theory  to 
complement  the  NTRM  corn  growth  model  for 
top  growth  and  total  dry  matter 
production.   The  unsaturated  flow  model 
(Dutt  et  al.  1972,  Shaffer  and  Gupta 
1981) ,  originally  based  on  works  reported 
by  Hanks  and  Bower  (1962) ,  was  updated  to 
include  the  layered  case,  and  upper 
boundary  conditions  were  included  that 
were  more  suitable  for  rapid  water  fluxes 
and  surface  evaporation.   The  evaporation 


and  infiltration  submodels  developed  by 
Linden  (1979)  were  incorporated  into  the 
overall  model. 

A  second  nitrogen  transformation  submodel 
developed  by  Molina  et  al .  (1982)  (see 
also  Chapter  7  of  this  handbook) ,  which 
is  more  detailed  and  less  empirical  than 
the  Shaffer  N- transformation  model  in  the 
areas  of  crop  residue  and  soil  organic 
matter  transformations  was  incorporated 
into  the  model.   This  submodel  also 
includes  the  options  of  using  Shaffer's 
transition  state  submodel  for 
nitrification  (Shaffer  et  al.  1977)  and 
an  expression  for  urea  hydrolysis 
previously  reported  by  Shaffer  et  al . 
(1969).   N-1-  -tracer  capabilities  were 
also  incorporated  into  the  new  N- 
transformation  model  and  the  overall  NTRM 
model . 

In  the  areas  of  crop  residue  and  tillage, 
new  submodels  and  information  transfers 
were  developed  to  simulate  the 
appropriate  interactions.   Capability  was 
added  to  allow  decay  of  residues  on  the 
soil  surface  and  interaction  of  these 
residues  with  crop  growth,  heat  flow,  and 
waterflow.   Similar  capabilities  were 
also  included  for  residues  incorporated 
into  the  soil. 

The  effects  of  tillage  practices  on  soil 
physical,  chemical,  and  biological 
properties  are  simulated  in  the  model. 
The  effects  of  those  changes  on  crop 
growth  are  then  modeled  using  a  separate 
set  of  relationships . 

The  soil  temperature  model  reported  by 
Gupta  et  al .  (1981)  was  incorporated  into 
the  NTRM  model  system.   This  submodel  is 
operated  as  an  independent  program  with 
input  to  the  interactive  model  via  data 
file  transfers. 

Capability  was  retained  in  the  overall 
NTRM  model  to  generate  soil  leachate 
volumes  and  constituent  concentrations 


suitable  for  input  to  aquifer  models  such 
as  the  one  reported  by  Ribbens  and 
Shaffer  (1976). 
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Chapter  2 

MODEL  STRUCTURE  AND  SUBMODEL  INTERACTIONS 

M.J.  Shaffer1 


The  basic  structure  of  the  NTRM  model  is 
shown  in  figure  2.1.   The  major 
computational  submodels  include  soil 
temperature,  unsaturated  flow, 
infiltration,  transpiration,  evaporation, 
tillage,  surface  residues,  nitrogen 
transformations,  chemical  equilibria, 
salt  transport,  crop  nutrient  uptake, 
crop  growth,  and  root  growth.   In 
addition,  there  are  two  executive 
subprograms  (EXECUTE  AND  COMBINE)  that 
control  the  logic,  flow,  and  integration 
of  the  computations ;  a  main  input 
subprogram;  and  a  number  of  smaller 
subroutines  (not  shown)  with  specific 
functions  including  data  output,  data 


transformations,  file  positioning,  rate 
calculations,  mass  balance  checking  and 
maintenance,  and  other  duties. 

The  general  input  and  output 
configuration  of  the  NTRM  model  is 
illustrated  in  figure  2.2.   A  detailed 
discussion  of  the  contents  of  each  file 
can  be  found  in  the  NTRM  User's  Guide 
(Shaffer  and  Pierce,  in  press).   Note 
that  the  soil  temperature  and  aquifer 
submodels  communicate  with  the  main 
interactive  soil  model  via  computer  files 
with  data  flow  in  only  one  direction. 
Substitution  of  alternative  submodels  for 
these  routines  can  be  done  by  providing 
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NTRM  model- -general  flowchart. 
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Figure  2.2. 

NTRM  model- -general  input  and  output 

configuration . 

programs  that  either  generate  or  accept 
the  appropriate  data  files.   Substitution 
within  the  NTRM  interactive  block  can  be 
done  with  general  input  control  switches 
or  by  replacing  submodels  or  equations. 
The  latter  case  requires  a  good  working 
knowledge  of  the  internal  source  code  and 
interfacing  conventions.   The  model  was 
written  in  a  modular  submodel  format  to 
make  this  task  easier.   Internal 
communication  is  via  common  blocks  and 
parameter  call  lists.   A  programmer's 
manual  will  be  published  to  assist  with 
interpretation  of  the  NTRM  FORTRAN  77 
source  coding. 

Output  is  generated  at  various  locations 
within  the  interactive  model  and  is 
written  to  disk  or  tape  files  for 
plotting  or  listing.   Leachate  volumes 
and  soluble  constituent  masses  are  output 
to  a  disk  or  tape  file.   These  data 
become  direct  input  to  the  saturated 
chemistry  submodel  (SATCHEM)  and  related 


submodels  that  simulate  chemical  and 
waterflow  processes  in  the  aquifer. 
Information  on  the  overall  aquifer  model 
(including  SATCHEM)  can  be  found  in 
Shaffer  et  al .  (1977).   Figures  2.3  and 
2.4  illustrate  the  general  inputs  and 
outputs  for  each  major  submodel. 

Processing  of  the  input  data  occurs 
within  each  routine,  which  then  generates 
new  information  or  modifies  existing 
data.   Details  of  the  internal  structures 
of  each  major  computational  submodel  are 
discussed  in  the  appropriate  chapters. 
Close  study  of  figures  2.1  through  2.4  is 
recommended  to  gain  a  good  appreciation 
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Figure  2.4. 

Submodel  usage  and/or  calculation  of 

chemical  and  biological  properties. 

for  the  general  information  flow 
throughout  the  NTRM  model  system. 

The  NTRM  model  is  written  in  FORTRAN  77 
computer  language  for  Control  Data 
Corporation  (CDC)  CYBER,  Cray  Research 
CRAY-1A,  IBM  mainframe,  IBM  PC/XT/AT,  and 
Hewlett-Packard  A700  machines.   A  version 
of  the  model  source  code  meeting  1977 
ANSI  FORTRAN  standards  is  available  to 
outside  users  on  request.   In  addition, 
an  extensive  operating  system  (Shaffer 
and  Pierce,  in  press)  for  the  NTRM  model 
has  been  developed  that  makes  it 
relatively  easy  to  access  and  use  the  CDC 
CYBER  and  other  versions  of  the  model. 
An  outside  user  can  enter  the  CDC  CYBER 
system  from  a  remote  time -share  terminal 
and  operate  the  model.   Results  can  be 
dumped  to  the   user's  terminal  or  the 
facilities  at  the  University  of 
Minnesota.   Separate  chapters  in  the 
user's  guide  (Shaffer  and  Pierce,  in 
press)  describe  this  system  in  detail. 


Storage  requirements  for  the  NTRM  model 
on  the  University  of  Minnesota  CYBER  74 
are  about  150K  binary  words  of  central 
memory.   The  model  uses  18  disk  or  tape 
files  during  execution.   Execution  times 
for  a  growing  season  vary  with  the  number 
and  types  of  submodels  being  used.   With 
the  entire  model  turned  on,  run  times  for 
a  growing  season  should  not  exceed  4  to  5 
seconds  (CRAY-1A  computer)  for  about 
$2.50.   Equivalant  run  times  on  the  CDC 
CYBER,  HP  A700,  and  IBM  AT  are  60 
seconds,  20  minutes,  and  30  minutes, 
respectively.   Long-term  and  other  runs 
using  special  reduced  configurations  of 
the  model  can  be  made  for  about  $0.50  or 
less  per  season  on  the  CRAY-1A. 

Some  common  model  configurations  that 
simplify  model  usage  are  given  in  Chapter 
4  "Model  Applications  and  Simplified 
Versions . " 


NTRM  Model  Spacial  Configuration 

The  horizon,  soil  segment,  and  nodal 
structures  of  the  NTRM  interactive 
submodel  block  (fig.  2.2)  are  illustrated 
in  figure  2.5.   The  model  has  a  one- 
dimensional  structure;  it  is  sectioned  by 
soil  depth  using  a  variable  spacing  for 
the  horizons,  segments  and  nodes. 
Horizons  contain  the  basic  soils  data  as 
input  into  the  model.   Segments  are  the 
computational  units  within  the  model. 
There  must  be  at  least  one  segment  for 
each  horizon.   Horizons  can  be  subdivided 
into  equally  spaced  segments  with  the 
same  DELX  (delta  X) .   When  a  simulated 
tillage  event  occurs  that  modifies  the 
soil  bulk  density,  the  model's  segment 
boundaries  are  modified  to  keep  the  soil 
mass  constant  within  each  segment.   The 
assumption  is  made  that  the  bottom  of  the 
soil  column  remains  fixed  while  other 
segment  and  horizon  boundaries  "float" 
with  their  positions,  depending  on  the 
transient  bulk  densities  of  each  segment. 
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Horizon,  segment,  and  nodal  structure 

within  the  one -dimensional  NTRM  model; 

equally  spaced  segments  and  nodes  in  each 

horizon. 


Nodes  are  the  computational  points  within 
the  unsaturated  flow  model.   Nodes  must 
be  placed  (among  other  locations)  on 
segment  and  horizon  boundaries . 
Individual  horizons  can  be  subdivided 
into  equally  spaced  nodes. 

A  surface  segment  or  compartment  is 
designed  to  accommodate  surface 
applications  of  crop  residues, 


fertilizers,  and  water.   Nitrogen  and 
carbon  transformations,  evaporation, 
leaching,  and  heat  transfer  are 
considered  in  conjunction  with  this 
surface  layer. 

When  tillage  is  expected  during  a 
simulation  run,  one  segment  boundary 
should  be  selected  to  correspond  to  the 
depth  of  the  tillage  operation. 

The  lower  boundary  corresponds  to  a  fixed 
water-table  position  or  to  a  point  in  the 
unsaturated  zone  with  a  relatively 
constant  water  content. 

The  nodal  structure  of  the  soil 
temperature  submodel  is  shown  in  figure 
2.6.   This  submodel  uses  a  finite 
difference  solution  to  the  heat  flow 
equation  written  in  one  dimension.   The 
nodal  configuration  is  similar  to  that  of 
the  unsaturated  flow  submodel,  except 
that  these  nodes  do  not  necessarily  have 
to  be  placed  on  soil  horizon  or  segment 
boundaries.   In  fact,  the  soil 
temperature  nodal  spacings  can  be 
independent  of  the  soil  horizons, 
segments ,  and  nodes  used  in  the 
interactive  subblock,  provided  the  soil 
temperature  nodes  extend  at  least  to  the 
bottom  boundary  of  the  unsaturated  zone 
as  defined  in  the  interactive  portion. 

Conversion  of  soil  temperatures  computed 
at  specified  nodes  to  soil  temperatures 
on  a  segment  basis  (finite  volume  of 
soil)  is  done  by  utility  routines  in  the 
interactive  submodel  block. 

The  two-dimensional  spacial  configuration 
of  the  aquifer  submodel  is  described  in 
the  appropriate  publications  on  this 
model,  for  example,  Ribbens  (1970). 


Submodel  Interactions 

With  the  exception  of  temperature 
calculations  and  the  aquifer  subprograms , 
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Nodal  structure  within  the  soil 

temperature  submodel. 

the  NTRM  submodels  are  entirely 
interactive,  with  numerous  feedback  loops 
operating  ( fig .  2.7). 

Top  and  Root  Growth 

Details  of  the  feedback  loops  involving 
top  and  root  growth  are  shown  in  figure 
2.7.   Note  that  a  minimum  number  of 
parameters  are  used  to  actually  operate 
the  feedback  systems.   The  basic 
parameters  involved  are  DEFNIT  for 
nitrogen,  RATI01  for  soil  water,  and 
RATIO  for  root  stress  effects.   Top 
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Figure   2.7. 

Interrelationships  of  top  and  root 

growth- -NTRM  model. 

growth  and  root  growth  are  considered 
simultaneously  in  the  model.   Nitrogen 
deficiencies  in  the  plant  that  lead  to 
reduced  growth  are  reflected  in  the 
DEFNIT  parameter.   Growth  is  reduced  in 
both  the  tops  and  roots  in  this  case. 
Water  deficiencies  in  the  root  zone  that 
lead  to  reduced  growth  are  reflected  in 
the  RATI01  parameter.   Again,  both  the 
tops  and  roots  show  reduced  activity 
addition,  other  conditions  in  the  root 
zone  that  limit  overall  root  and  top 
development  (any  possible  root 
distribution  changes  taken  into  account) , 
such  as  soil  temperature,  strength, 
salinity,  aeration,  nitrogen,  and 
suction,  and  that  are  not  already 
accounted  for  in  the  DEFNIT  and  RATI01 
parameters  are  included  in  the  RATIO 
parameter.   These  can  be  thought  of  as 
indicators  of  situations  where  specific 
root  zone  conditions  will  not  allow  the 
development  or  redistribution  of  the 
overall  root  mass  as  projected. 


Definitions  of  RATI01 
Parameters 


DEFNIT,  and  RATIO 


RATI01  is  defined  as  the  ratio  of  actual 
to  potential  water  extraction  by  the 
roots  over  a  time  step.   Here,  potential 
water  extraction  takes  into  account  the 
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physiological  state  of  the  crop  and  stage 
of  development .   Only  if  the  crop  runs 
out  of  water  on  this  basis  will  RATI01  be 
less  than  1.0.   Permanent  wilt  is  taken 
as  the  lower  limit  for  water  extraction. 
Potential  water  extraction  as  a  function 
of  depth  is  assumed  to  be  proportional  to 
current  root  distribution. 

DEFNIT  is  defined  as  the  ratio  of  actual 
to  potential  uptake  of  NO3" -N+NH4+-N  by 
the  crop  for  that  time  step.   N03_-N  and 
soluble  NH^+-N  are  treated  as  being 
equally  accessible  by  the  roots. 
Potential  uptake  of  N  is  based  on  the 
current  physiological  state  of  the  plant 
and  the  growth  stage.   Potential 
extraction  of  N  as  a  function  of  depth  is 
assumed  to  be  proportional  to  actual 
water  uptake  with  depth.   When  a  nitrogen 
stress  condition  is  removed,  a  lag  period 
of  up  to  a  few  days  is  used  before  full 
nitrogen  uptake  resumes. 

RATIO  is  defined  as  the  ratio  of  actual 
to  potential  production  of  root  mass  for 
a  time  step.   Potential  root  mass  is 
defined  as  the  mass  of  roots  that  the 
plant  could  produce  based  on  its 
physiological  state  and  its  growth  stage. 
This  includes  prior  consideration  of  the 
DEFNIT  and  RATI01  parameters. 


Tillage  and  Soil  Property  Interactions 

Interactions  between  tillage  and  soil 
properties  are  illustrated  in  figure  2.8. 
Note  that  transient  values  for  soil 
properties  affect  the  impact  tillage  has 
on  the  soil  properties.   The  model  either 
modifies  the  soil  properties  according  to 
input  specifications,  or  it  uses  internal 
submodels  to  estimate  the  changes  that 
would  occur.   This  applies  to  individual 
tillage  events  as  well  as  to  more  gradual 
changes  that  occur  between  events. 

More  specific  data  transfers  between  the 
tillage  submodel  and  other  routines  are 
shown  in  figure  2.9.   Note  that  tillage 
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Figure  2.8. 

Interactions  of  tillage  with  soil 

properties . 

events  (and  between- tillage  events) 
directly  affect  several  submodels.   Other 
submodels  are  affected  indirectly  because 
of  the  changing  of  related  processes. 
The  TILL  parameter  for  increased 
biological  activity  and  updated  values 
for  soil  nitrogen  concentrations  and  soil 
water  content  are  provided  for  the  carbon 
and  nitrogen  transformation  submodel. 
Soil  water  contents  and  unsaturated  flow 
properties  are  updated  for  use  by  the 
unsaturated  flow  submodel.   Values  for 
soil  bulk  density,  soil  water  content, 
water  content  at  saturation  (0S) ,  matric 
potential,  osmotic  potential,  soil 
nitrogen  concentrations,  and  root  pruning 
are  recomputed  for  input  into  the  root 
growth  submodel.   Residue  mass,  percent 
cover,  surface  reflection  coefficients, 
and  surface  roughness  are  updated  for  the 
evaporation  submodel.   Values  for  soluble 
and  exchangeable  constituent 
concentrations  are  modified  for  use  in 
the  chemical  equilibria  submodel. 


Soil  Nitrogen  Interactions 

Interaction  of  the  nitrogen  and  carbon 
transformation  submodel  with  other 
submodels  is  shown  in  figure  2.10.   The 
COMBINE  executive  submodel  performs  the 
numerical  integration  of  the  overall 
system  with  time.   Note  that  crop  and 
root  growth,  root  slough,  nitrogen  and 
carbon  transformations,  ion  exchange 
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Figure  2.9. 

Exchange  of  tillage  data  between 

submodels . 

involving  NH^-N,  and  solute  transport 
are  all  considered  simultaneously. 
Constituent  masses  for  the  nitrogenous 
species  as  a  function  of  depth  are  stored 
in  a  common  program  block  area  in  the 
model.   This  allows  access  by  all  the 
submodels  as  required.   For  a  given  time 
step,  each  submodel  uses  the  same 
nitrogen  mass  data.   Updates  are  made 
after  each  submodel  has  completed  its 
calculations.   Any  mass  overruns 
(negative  values)  are  redistributed  back 
over  the  processes  involved  in  order  to 
maintain  mass  balance  in  the  system. 

The  solute  transport  submodel  uses 
soluble  NH^+-N  and  N03~-N  masses,  soil 
water  contents,  and  water  fluxes  as  a 
function  of  depth  to  compute  transport  of 
these  species  (mass/segment/time  step) . 
The  chemical  equilibria  submodel  uses 
current  values  for  soluble  and 
exchangeable  NH4  -N  concentrations  along 


with  similar  values  for  Ca"1-*",  Na+,  and 
Mg"1-1",  and  concentrations  of  HCO3",  CC>3=, 
SC>4=,  CI",  and  N03~-N  to  compute  updated 
values  for  exchangeable  and  soluble 
NH4+-N.   The  root  growth  submodel  uses 
NH4+-N  and  N03~-N  concentrations  and  the 
crop  nitrogen  deficiency  parameter 
(DEFNIT)  in  its  calculations  and  returns 
the  daily  root  slough  for  recycling  by 
the  nitrogen  and  carbon  transformation 
submodel.   The  crop  growth  submodel  uses 
the  DEFNIT  parameter  and  other  variables 
to  compute  potential  uptake  of  soil  N. 
NH4+-N  and  NC>3"-N  concentrations  and 
other  information  are  used  to  compute 
actual  N  uptake. 


Soil  Water  Interactions 

Submodel  interactions  involving  soil 
water  are  illustrated  in  figure  2.11. 
The  unsaturated  flow  submodel  uses  soil 
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water  contents  and  flow  properties  to 
compute  updated  values  for  soil  water 
contents,  matric  potentials,  and  fluxes. 
Soil  water  contents  are  used  in  the  root 
growth,  carbon  and  nitrogen 
transformation,  chemical  equilibria, 
solute  transport,  and  tillage  submodels. 
Soil  suctions  are  used  in  the  root  growth 
submodel,  while  flow  properties  are  input 
to  the  tillage  submodel,  and  soil  water 
fluxes  are  used  by  the  solute  transport 
submodel.   Submodel  feedback  based  on  the 
results  of  using  the  soil  water 
information  takes  the  form  of  crop  and 
root  masses  produced,  transformation 
rates  for  carbon  and  nitrogen,  movement 
of  soluble  constituents,  the  soil 
chemical  equilibrium  status,  and 
potential  values  for  evaporation  and 
transpiration.   The  tillage  submodel 
makes  updates  in  soil  properties  based, 
in  part,  on  information  provided  from  the 
unsaturated  flow  submodel. 
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Chapter  3 

OVERALL  MODEL  VALIDATION  AND  VERIFICATION 

M.J.  Shaffer1 


Development  of  the  NTRM  model  over  time 
is  illustrated  in  figure  3.1.   Note  that 
progress  on  the  model  was  rapid  during 
1980  and  1981  with  a  base  having  been 
previously  established  during  the  1960 's 
and  1970' s.   This  progress  included  both 
validation  and  verification,  activities 
on  the  individual  submodels  (see 
appropriate  submodel  chapters)  and  the 
overall  integrated  model. 

With  complex  models  such  as  NTRM, 
absolute  validation  and  verification  is 
generally  only  a  theoretical  possibility 
that  occurs  at  infinite  time  in  the 
future.   There  is  always  a  finite 
probability  that  the  model  still  contains 
bugs  or  that  the  underlying  concepts  are 
inadequate  or  inappropriate. 

The  purpose  of  continued  validation  and 
verification  is  to  minimize  the  risks  to 
the  extent  practical.   This  usually  means 
(based  on  experience  with  other  models) 
that  the  model  developers  will  debug  and 
verify  their  models  to  the  limits  of 
their  budgets  or  interest  and  that  model 
users  will  identify  and  correct  (often 
assisted  by  the  developers)  any  remaining 
basic  model  deficiencies.   Experience  has 
shown  that  user  feedback  is  one  of  the 
most  effective  and  efficient  ways  to  make 


progress  with  model  validation  and 
verification  once  a  reasonable  level  of 
confidence  has  been  reached. 

The  present  state  of  the  NTRM  model  is 
about  that  shown  in  figure  3.1.   Some 
additional  validation  and  verification 
will  be  done  by  the  NTRM  development 
group.   The  user  community  will  furnish 
the  rest.   A  staff  will  be  maintained  to 
assist  the  users  with  any  model  problems 
that  may  develop.   The  details  of  the 
model  validation  and  verification 
completed  for  the  overall  model  and  some 
submodels  are  available  in  various  NTRM 
related  publications. 

The  following  summary  is  presented  to 
give  the  reader  an  appreciation  for 
overall  model  validation  and  verification 
efforts.   Model  validation  in  terms  of 
code  checking  and  testing  and  internal 
mass  balance  checks  has  been  done  at  each 
stage  of  submodel  development  and  after 
incorporation  of  each  submodel  into  the 
overall  NTRM  system.   Members  of  the 
research  team  have  examined  in  detail 
portions  of  the  model  developed  by  other 
members .   Errors  have  been  corrected  as 
soon  as  they  were  identified.   Model 
sensitivity  runs  have  been  done 
extensively  to  determine  if  simulated 
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Figure  3.1. 

NTRM  model  development  and  verification. 


1  Research  soil  scientist,  U.S.  Department 
of  Agriculture,  Agricultural  Research 
Service,  Soil  and  Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 
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output  appears  reasonable  over  the 
expected  range  of  input  parameters. 
Graphical  output  has  been  used  to  a  real 
advantage  in  these  analyses. 

Field  testing  of  the  overall  integrated 
model  has  been  a  continuing  process  as 
additional  data  sets  are  analyzed.   For 
example,  corn  yield  results  for  four 
different  locations ,  three  in  Minnesota 
and  one  in  Ohio,  are  presented  in  figure 
3.2.   Yield  differences  seen  here  are 
primarily  the  result  of  climatic 
differences  as  opposed  to  variations  in 
management . 

Figures  3.3  and  3.4  are  examples  of  crop 
dry  matter  production  and  crop  nitrogen 
uptake  curves ,  in  this  case  for  predicted 
and  observed  corn  data  collected  on  a 
control  site  at  Apple  Valley,  MN.   The 
model  run  indicated  that  nitrogen  was  not 
limiting  and  a  relatively  shallow  root 
system  would  be  expected  because  of  ample 
nutrient  availability  and  moist  soil 
conditions  at  shallow  depths.   This  was 
confirmed  by  root  samples  in  the  field. 
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Figure   3.3. 

Crop  dry  matter  production  for  predicted 
and  observed  corn  data  collected  on  a 
control  site  at  Apple  Valley,  MN. 
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NTRM  model- -corn  yields. 
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Figure  3.4. 

Crop  nitrogen  uptake  for  predicted  and 
observed  corn  data  collected  on  a  control 
site  at  Apple  Valley,  MN. 
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Chapter  4 

MODEL  APPLICATIONS  AND  SIMPLIFIED 

VERSIONS 

M.J.  Shaffer1 


Introduction 

The  NTRM  model  is  a  comprehensive 
simulator  of  the  crop-soil-water 
continuum.   The  model  contains 
considerable  detail  about  crop  growth  and 
yield  response  to  many  climatic  and  soil 
variables.   It  emphasizes  nitrogen, 
tillage,  and  residue  aspects  of  the 
overall  system. 

Application  of  this  tool  and  of 
information  derived  from  its  application 
is  possible  in  many  areas,  including 
research,  engineering,  forecasting, 
teaching,  extension,  and  direct  use  by 
farmers.   A  computer  operating  system  has 
been  developed  for  the  model  to 
facilitate  direct  operation  by  extension 
personnel,  researchers,  engineers,  and 
forecasters . 


Extension  and  Farm  Applications 

Extension  personnel  involved  in  the 
development  of  guidelines  and 
recommendations  will  find  the  model 
system  extremely  useful.   Model 
simulation  runs  can  be  made  in  lieu  of 
lengthy  and  costly  field  plot 
experiments  and  thus  greatly  improve  the 
effective  rate  at  which  problems  are 
solved.   The  model  would  be  extremely 
useful  in  cases  where  unusual  or  unique 
management  problems  occur  and  immediate 
answers  are  needed. 

Management  of  agriculture  on  a  shorter 
term  basis,  such  as  seasonal  or  monthly 
time  increments,  is  an  area  where  the 
NTRM  model  can  provide  answers .   The 
present  status  of  the  soil-water-crop 
system  can  be  assessed  in  terms  of 

1  Research  soil  scientist, 
U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 


available  management  techniques  and 
forecasted  environmental  conditions  to 
project  the  best  short-term  management 
scheme . 

The  model  can  be  used  to  directly  or 
indirectly  generate  information  in 
condensed  forms  such  as  simplified 
equations,  nomographs,  charts,  and 
written  text  for  use  in  farm  and 
extension  programs.   This  would  include 
use  in  computerized  extension  systems, 
such  as  Green  Thumb  and  Agnet,  and  in  the 
U.S.  Department  of  Agriculture, 
Agricultural  Research  Service  COFARM 
(Coordinated  Farm  and  Research  Manage- 
ment) model. 

For  agricultural  planning,  the  NTRM  model 
has  several  important  applications.   It 
enables  planners  to  assess  long-term 
impacts  of  physical,  chemical,  and 
biological  conditions  (whether  natural  or 
induced  by  human  activity)  on  crop 
yields,  soil  conditions,  aquifer  quality, 
and  surface  stream  quality.   These 
results  can  provide  important  input  into 
long-range  planning  decisions  involving 
agriculture.   The  model  can  also  enable 
planners  to  assess  the  long-term 
environmental  impact  of  irrigation 
projects,  conservation  practices  (for 
example,  no  tillage  and  conservation 
tillage),  energy  proposals  (for  example, 
crop  residue  usage  and  coal  and  oil  shale 
production) ,  and  other  management 
practices ,  such  as  usage  of  chemical 
fertilizers  and  sewage  sludge  and 
effluent . 

Research  Applications 

The  model  in  its  present  form  has  many 
research  applications  for  both  laboratory 
and  field  studies  where  extensive  data 
collection  is  cumbersome  or  impossible  or 
where  new  ideas  or  research  directions 
are  needed.   The  model  allows  rapid 
interpolation  to  generate  missing  data 
sets  in  large-scale  experiments.   In 
these  cases,  the  savings  in  time  and 
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money  can  be  significant.   And,  a 
particular  system  under  investigation  can 
be  rapidly  varied  or  manipulated  around  a 
particular  problem,  either  to  suggest  new 
research  or  to  predict  the  most 
profitable  lines  of  research. 

Analysis  of  model  sensitivity  has  a  broad 
range  of  research  applications.   Input 
parameters  and  model  coefficients  can  be 
varied  systematically  to  determine  which 
variables  have  the  most  significant 
impact  on  the  desired  output.   Planners 
can  then  design  a  data-collection  program 
that  emphasizes  the  most  significant 
input  variables  and  coefficients- -rather 
than  using  a  hit-and-miss  procedure  based 
on  someone's  experience  or  a  blanket 
approach  in  which  a  wide  spectrum  of  data 
is  collected. 

Sensitivity  analyses  can  help  indicate 
the  degree  of  detail  needed  to  accomplish 
the  objectives  in  a  particular  type  of 
simulation  analysis  or  overall  study. 
One  method  that  can  be  used  to  develop  a 
model  or  approach  that  falls  near  the 
optimal  is  to  simplify  an  overly  detailed 
model  by  using  sensitivity  analyses. 
Indeed,  the  approaches  being  used  in  most 
research  studies  could  benefit  from  model 
sensitivity  analyses. 

Teaching  Applications 

The  NTRM  model  is  an  excellent 
instructional  tool  in  the  area  of  the 
soil-water-crop  system.   Each  submodel 
can  be  taught  as  a  separate  unit ,  and  the 
interactions  between  the  submodels 
(subsystems)  can  be  examined  in  some 
detail.   Much  can  be  learned  about  the 
real  world  of  soil,  water,  and  crop 
interrelationships  by  manipulating  the 
input  data  and  observing  model  responses . 

The  best  way  to  illustrate  the 
capabilities  of  the  NTRM  model  in  the 
extension  area  is  by  examining  several 
examples.   In  most  cases,  the  model  would 
be  applied  where  there  is  a  deficiency  in 


existing  data  concerning  a  particular 
crop  response  or  other  problem.   An 
example  of  this  would  be  a  need  by 
farmers  to  know  the  yield  response  of 
corn  to  various  tillage  practices  for 
soils  and  climatic  conditions  not 
extensively  studied  at  the  experiment 
stations . 

Tillage 

A  typical  yield  response  nomograph  that 
can  be  generated  by  the  NTRM  model  is 
shown  in  figure  4.1.   The  input  data 
needed  by  the  model  to  produce  this 
output  are  outlined  in  the  tabulation 
below. 

NTRM  model  input  data  needed  for 
nomograph  generation2 

Climatic  data: 

1.  Daily  max-min  air  temperature 

2 .  Daily  wind 

3.  Daily  precipitation 

Soil  properties: 

1.  Bulk  density  (initial) 

2.  Soil  water  content  (initial) 

3.  Soil  water  characteristic  curves 

(optional) 

4.  Surface  roughness  and  reflection 

coefficients  (initial) 

5.  Soil  extract  analysis- -nitrogenous 

species  and  major  ions 
(optional) 

6.  Soil  temperatures  (initial) 

7 .  Cation  exchange  capacity 

(optional) 

8.  Soil  organic  nitrogen  and  C/N 

ratio 

9.  Soil  texture 

10.   Maximum  rate  of  infiltration 

2  Note:   This  list  is  based  on  a  number 
of  variables  being  set  at  default 
values  in  the  model.   These  can  be 
modified  by  the  user  for  research 
applications . 
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Management  practices 


100 


Crop  variety  (coefficients) 
Planting  date 

Tillage  dates  and  practice 
Nitrogen  fertilizer  and  crop 

residue  application  dates, 

amounts,  types,  depths,  and 

percent  of  coverage 
Irrigation  dates  and  amounts  (if 

any) 
Salinity  of  irrigation  water  (if 

any) 
Starting  date  of  run 


If  the  graph  in  figure  4.1  were  based  on 
10  years  of  weather  records  and  simulated 
three  tillage  practices,  30  computer  runs 
would  be  needed,  at  a  total  cost  of  about 
$75.   Use  of  an  average  or  specific  set 
of  weather  data  would  reduce  this  to 
$7.50.   Optional  application  of  less 
scientifically  rigorous  submodels  within 
the  NTRM  system  could  produce  the 
nomograph  for  as  little  as  $1.00.   Once 
the  input  data  are  assembled,  the 
computer  runs  can  be  made  in  a  single 
night.   These  dollar  and  time  costs 
compare  very  favorably  with  the  10 -year 
period  and  many  years  of  effort  needed  to 
duplicate  this  with  field-plot  research. 

In  addition,  once  the  nomograph  is 
available,  extension  agents  can  make  it 
available  to  many  different  farmers  with 
similar  conditions  by  placing  it  on  an 
extension  computer  system.   The  cost  to 
each  farmer  for  this  information  would 
amount  to  the  computer  connect  charges 
and  extremely  small  central  processing 
costs  needed  to  list  the  graphic  on  the 
farmer's  own  terminal  or  an  extension 
office  terminal. 

Nitrogen 

For  nitrogen  fertilizer  recommendations, 
nomographs  derived  from  the  NTRM  model 
can  be  applied  extensively.   An  example 
of  fertilizer-related  output  is  shown  in 
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Figure  4.1. 

Corn  yield  response  to  climate  and 

tillage. 

figures  4.2-4.5.   Here  a  fertilizer 
recommendation  can  be  derived  as  a 
function  of  a  number  of  different 
parameter  values  specific  for  each  soil. 
By  knowing  the  nitrogen  requirements  of 
the  crop  and  the  initial  soil  nitrogen, 
the  user  can  apply  these  four  graphs  in 
combination  to  arrive  at  an  estimate  of 
the  fertilizer  required  for  an  individual 
soil.   In  special  cases,  individual 
nomographs  can  be  used  to  decide  when  any 
additional  fertilizer  may  be  needed. 
These  nomograph  sets  can  be  custom 
tailored  for  individual  crops,  soils,  and 
locations . 

In  addition  to  fertilizer 
recommendations,  nomographs  relating  to 
environmental  aspects  of  soil  nitrogen 
(for  example,  N03~-N  leaching)  can  be 
used  for  studies  on  pollution  abatement 
and  environmental  impact. 
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Figure  4.3. 

Denitrif ication  after  a  6-cm  rain. 
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Figure  4.4. 

Soil  N  mineralization   in  a  season. 
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Figure  4.5  illustrates  a  nomograph  that 
gives  the  additional  fertilizer  required 
(or  added)  for  various  crop  residues 
incorporated  into  the  soil .   Other 
nomographs  are  possible  based  on  output 
from  the  NTRM  model.   For  example,  a 
graphic  relating  changes  in  soil  surface 
temperatures  for  different  crop  residues 
could  be  prepared.   Similar  nomographs 
could  also  be  constructed  for 
infiltration,  residue  decay  products,  and 
soil  erosion. 

From  the  above  examples ,  it  is  evident 
that  a  wide  range  of  useful  nomographs  is 
possible  based  on  output  from  the  NTRM 
model.   Once  management  problems  are 
defined  in  specific  areas,  suitable 
computer  runs  can  be  designed  to  provide 
nomographs  or  other  output  to  aid  in 
solving  these  problems  on  a  practical 
level. 

Introduction  of  these  aids  into  an 
interactive  computer  network  will  make 
this  information  available  to  a  wide 
variety  of  users.   A  sample  nomograph 
listed  on  a  Texas  Instruments  745 
portable  terminal  is  illustrated  in 
figure  4.6.   This  graphic  was  produced  in 
less  than  20  seconds  at  a  cost  of  a  few 
cents.   The  graphic  can  be  listed  on  any 
time -share  terminal  with  at  least  an  80 
column  line  printer.   This  illustrates 
the  relative  ease  and  efficiency  by  which 
information  generated  via  a  highly 
sophisticated  scientific  tool  can  be  made 
available  for  general  use.   This  short 
transfer  time  greatly  reduces  the  lag 
between  research  findings  and  their 
practical  application. 


The  NTRM  model  contains  a  number  of 
control  switches  that  allow  the  user  to 
select  submodel  configurations  that  are 
appropriate  for  a  particular  problem.   It 
may  not  be  necessary  or  desirable  to 
operate  all  the  NTRM  model  capabilities 
during  a  particular  simulation. 

Table  4.1  is  a  list  of  some  common 
simulations  that  the  NTRM  model  can 
handle  without  modification.   One 
important  type  of  application  for  certain 
configurations  is  the  ability  to  conduct 
sensitivity  analyses  that  compare  results 
obtained  for  various  conditions  with 
results  for  the  ideal  case.   This  allows 
a  rapid  assessment  of  a  particular 
treatment  with  respect  to  its  maximum 
possible  benefits.   Different  processes 
can  be  compared  to  see  which  ones  would 
or  might  have  the  greatest  impact  on  the 
results . 

It  should  be  possible  to  identify  the 
most  efficient  model  configuration  for  a 
given  application.   What  is  needed  is  the 
simplest  model  design  that  will  give  the 
required  accuracy.   It  is  relatively  easy 
to  operate  a  more  complex  model  like  NTRM 
in  a  variety  of  configurations  to  help 
identify  this  design.   However,  a  model 
that  is  extremely  simple  to  start  with 
cannot  be  tested  beyond  its  capabilities 
without  additions. 
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YIELD  RESPONSE  AS  A  FUNCTION  OF  CLIMATE  AND  TILLAGE 
BASE  YIELD  =  150  BU/AC 
BASE  TEMP  =  21°C 
BASE  WATER  CONTENT  =  SEE  WC  NOMOGRAPH 
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Figure  4.6. 

Sample  nomograph- -NTRM  model,  June  1981. 
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Table  4.1. 

Some  common  simplified  NTRM  model  configurations 


Configuration 


Description 


Chemical  Equilibria  Bypass 


The  chemical  equilibria  submodel 
is  not  called.   Major  cations  and 
anions  (other  than  nitrogenous 
species)  are  treated  as 
conservative . 


Ideal  Crop  Response  to  Nitrogen 


Ideal  Root  Response 


Ideal  Crop  Response  to  Water 


Nitrogen  Transformation  Bypass 


Steady- State  Waterflow 


The  crop  is  grown  with  a 
nonlimiting  response  to  nitrogen 
deficiencies  even  though  nitrogen 
may  be  depleted  in  the  root  zone . 

The  crop  roots  are  grown  as  though 
ideal  conditions  exist  for  root 
development. 

The  crop  is  grown  as  though  supply 
is  nonlimiting. 

Nitrogenous  species  are  treated  as 
being  conservative  with  respect  to 
biological  N  transformations. 

Water  movement  within  the  soil 
profile  occurs  at  a  constant  rate 
and  water  content. 


Simplified  Potential  Evaporation 


Soil  Temperature  Data  Read 


Tillage  Data  Read 


Potential  evaporation  from  the 
soil  surface  is  computed  using  a 
simple  relationship  between  LAI 
(leaf  area  index)  and  pan 
evaporation. 

Soil  temperature  data  are  read 
weekly. 

Information  on  changes  in  the  soil 
due  to  tillage  is  read  in  from 
cards . 
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Chapter  5 

FUTURE  RESEARCH  NEEDS 

M.J.  Shaffer1 


Large  model  systems  similar  to  NTRM  will 
be  refined  and  expanded  as  new  knowledge, 
needs ,  and  computer  hardware  are 
developed.   Current  model  designs  should 
be  thought  of  as  plateaus  that  serve 
their  purpose  until  new  models  become 
available . 

Several  parts  of  the  NTRM  model  need 
future  development.   We  have  research 
scheduled  to  address  these  areas .   The 
chemical  equilibria  submodel  will  be 
expanded  to  include  acid  soils.   The 
submodel  is  presently  limited  to  alkaline 
and  neutral  soil  conditions.   Submodels 
will  be  developed  or  adapted  for 
processes  associated  with  phosphorus, 
potassium,  and  trace  metals.   This  will 
enlarge  the  scope  of  NTRM  management 
capabilities.   Additional  crops 
(soybeans,  small  grains,  pasture,  and 
others)  are  being  added.   We  are 
consolidating  crop  growth  processes  into 
common  subroutines  to  avoid  duplication 
between  crop  types.   Submodels  will  be 
added  that  simulate  the  effects  of  crop 
diseases  and  pests.   These  routines  will 
be  made  interactive  with  the  rest  of  the 
model.   A  submodel  will  be  added  to 
simulate  exchange  of  gases  between  the 
soil  and  the  atmosphere.   This  will 
improve  predictive  capabilities  in 
several  submodels.   Model  validation 
and  verification  will  continue  with 
respect  to  the  overall  model  and  the 
individual  submodels.   The  geographical 
range  of  data  used  in  these  efforts  will 
be  expanded  outward  from  Minnesota  and 
the  ARS  North  Central  Region.   The 
results  will  be  published  in  appropriate 
professional  journals  and  reports.   It  is 
expected  that  considerable  effort  will  be 
expended  in  this  area. 

1  Research  soil  scientist, 
U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 


Stochastic  and  geostatistical  methods  of 
data  analysis  will  be  adapted  or 
developed  for  use  with  NTRM  model  input 
data  and  coefficients.   This  will  allow 
the  assignment  of  well-defined  confidence 
limits  to  NTRM  predicted  output.   The 
area  of  management  response  forecasting 
should  benefit  from  this  work. 

Reduced  or  simplified  versions  of  the 
NTRM  model  will  be  developed  to  operate 
on  the  new  generation  of  microcomputers 
such  as  the  HP  A700,  IBM  PC/XT/AT,  APPLE 
Macintosh,  TERAK,  and  RADIO  SHACK.   We 
will  try  to  retain  the  maximum  possible 
NTRM  accuracy  within  the  hardware  limits 
of  these  machines.   This  should  greatly 
expand  the  public  availability  of  NTRM 
software . 

The  NTRM  model  will  be  used  to  make 
regular  forecasts  in  the  area  of 
nitrogen- tillage -residue  management. 
These  forecasts  in  the  form  of  maps, 
nomographs ,  graphs ,  and  text  will  be  made 
available  to  the  public  through 
conventional  media  or  on  the  computer 
information  network  systems  such  as  AGNET 
and  GREEN  THUMB.   An  interactive  user 
system  for  nitrogen,  tillage,  and  residue 
management  (COFARM)  has  been  developed 
at  St.  Paul.   This  system  uses  output 
from  the  NTRM  model  along  with  existing 
information  and  data  bases ,  including 
soils  and  climate,  to  provide  the  farmer 
with  site -specific  management 
recommendations.   This  system  is  being 
field- tested  in  southeastern  Minnesota. 
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Chapter  6 

SOIL  TEMPERATURE  SUBMODEL 

S.C.  Gupta,  J.K.  Radke,  and  W.E.  Larson1 


Symbols 

The  symbols  used  most  often  in  this 
chapter  are  defined  below.   Others  are 
defined  in  the  text. 

A-p  =  difference  in  amplitude  of  the 

sine  function  between  soil 

surface  and  air  temperature 
D  =  difference  between  predicted 

and  measured  temperatures 

averaged  over  the  growing 

season,  °C 
N  =  number  of  observations 
P  =  period,  hours 
SD  =  standard  deviation  of  D,  °C 
T  =  soil  temperature,  °C 
TQ  t  =  upper  boundary  temperature , 

°C 
Tl  t  =  lower  boundary  temperature , 

°C 
Tj4  =  measured  temperature,  °C 
Tp  =  predicted  temperature,  °C 

hourly  air  temperature  at  2-m 

height,  °C 

minimum  temperature  of  air  at 

2-m  height,  °C 

maximum  temperature  of  air  at 

2-m  height,  °C 

minimum  temperature  under  the 

soil  surface,  °C 

maximum  temperature  under  the 

soil  surface,  °C 

normalized  hourly  air 

temperature 

averaged  normalized  hourly  air 

temperature 

1  Gupta  is  an  associate 
professor,  Soil  Science 
Department,  University  of 
Minnesota,  St.  Paul,  MN  55108; 
Radke  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Rodale 
Research  Center,  Kutztown,  PA 
19530;  and  Larson  is  head, 
Soil  Science  Department, 
University  of  Minnesota,  St. 
Paul,  MN  55108. 


La,  t 

T   = 
■"-an 

T   = 
^ax 

T   = 
^sn 

T   = 

isx 


Ta,  t 

T? 


t  p  =  standard  deviation  of  the 

averaged  normalized  hourly  air 
temperature 
UBT  =  upper  boundary  temperature , 
°C 

Z  =  soil  depth,  cm 

i  =  depth  increment,  cm 

j  =  time  increment,  hours 

p  =  probability  level 

t  =  time ,  hours  or  days 
t]_,t2,t3  =  time  (hours  when  the  air 

temperature  is  equal  to  or 

slightly  higher  than  the 

reference  temperature 

thermal  diffusivity  cm^/h 

t-t2 

growing  degree  days 


a  = 

t„  = 


'X 

GDD 


Introduction 

Reduced  tillage  or  no- till  systems 
commonly  leave  all  or  part  of  the 
residues  from  the  previous  crop  on  the 
soil  surface.   These  residues  are 
important  in  helping  control  erosion  and 
evaporation.   Surface  residues  left  after 
a  tillage  operation  influence  the  heat 
and  water  balances  of  soil  profiles  over 
a  growing  season.   In  some  cases,  this 
influence  of  crop  residues  could  be 
critical  for  seed  germination,  plant 
growth,  and  crop  yield. 

The  effect  of  crop  residues  on  heat 
balance  of  soil  profiles  depends  largely 
on  the  quantity,  type,  and  placement  of 
these  residues.   Generally,  surface 
residues  depress  (Van  Wijk  et  al .  1959) 
the  soil  temperatures  at  the  seed  zone 
depth  by  (a)  acting  as  an  insulating 
layer  and  (b)  reflecting  a  greater 
fraction  of  solar  radiation  back  to  the 
atmosphere .   The  reflection  is  due  to  the 
lighter  color  of  crop  residues  compared 
with  that  of  the  soil  surface. 

Depressed  seed  zone  temperatures  have 
often  been  related  to  a  delay  in  seed 
germination,  slow  corn  seedling  growth, 
and  reduced  yield  of  corn  in  areas  with 
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cool  and  wet  springs  like  the  northern 
Corn  Belt  (Van  Wijk  et  al .  1959,  Burrows 
and  Larson  1962,  Olson  and  Schoeberl 

1970,  Griffith  et  al .  1973,  Olson  and 
Horton  1975) .   A  review  by  Willis  and 
Amemiya  (1973)  details  the  effects  of 
tillage  on  soil  temperature  and  the 
interaction  of  tillage  and  soil 
temperature  with  growth  of  various  crops . 
Splinter  (1974),  Duncan  (1975),  and 
Jenkins  (1975)  described  the  effects  of 
temperature  on  yield  of  various  crops 
through  computer  modeling  of  root  and  top 
growth.   The  NTRM  model,  like  some  other 
crop  growth  models,  needs  hourly  air 
temperature  as  input  for  calculation  of 
crop  growth  and  daily  air  temperature  as 
input  for  calculation  of  nitrogen  trans- 
formations and  root  growth.   The  need  for 
a  temperature  submodel  that  predicts  soil 
temperatures  from  parameters  (like  hourly 
air  temperature)  routinely  required  in 
crop  growth  models  necessitated  the 
undertaking  of  this  study. 

Two  types  of  soil  temperature  models  have 
been  discussed  in  the  literature: 
empirical  models  and  physical  models. 
Empirical  soil  temperature  models 
(Hasfurther  and  Burman  1974,  Toy  et  al . 
1978,  Cruse  et  al .  1980)  are  based  on  the 
statistical  relationships  between  soil 
temperature  at  some  depth  and 
climatological  and  soil  variables.   These 
models  are  easy  to  use,  but  the 
requirements  for  a  large  data  base  to 
develop  empirical  coefficients  limit 
their  application  to  treatments  and 
locations  other  than  those  used  in 
developing  the  coefficients. 

Physical  soil  temperature  models 
(Wierenga  and  deWit  1970,  Hanks  et  al . 

1971,  Gupta  et  al .  1981,  1982),  on  the 
other  hand,  are  based  on  the  principles 
of  heat  flow  in  soils  and  are  not  site  or 
treatment  specific.   However,  in  these 
models ,  knowledge  of  upper  boundary 
temperatures  is  necessary  to  predict  soil 
temperatures  with  depth  and  time. 


The  objectives  of  this  study  were  to  (1) 
develop  procedures  for  estimating  upper 
boundary  temperatures  for  bare  and 
residue  covered  soils  with  and  without  a 
growing  corn  crop  from  climatic  and  soil 
parameters  that  are  routinely  used  in 
crop  yield  prediction  models  and  (2) 
evaluate  the  usefulness  of  a  physical 
model  for  predicting  root  zone 
temperatures  when  upper  boundary 
temperatures  are  estimated  rather  than 
measured. 

Model  Description 

The  root  zone  temperature  model  is  based 
on  the  partial  differential  equation 
describing  heat  flow  in  soils. 


dT/dt=3(a3T/aZ)/dZ, 


[i: 


Equation  1  can  be  numerically  solved  for 
soil  temperature,  T,  with  depth,  Z,  and 
time,  t,  if  the  thermal  diffusivity,  a, 
along  with  the  initial 

T=T(Z)  [2] 

and  boundary  conditions 

Tot=T(t)  and  [3] 

TL|t=T(t)  [4] 

are  known. 

The  numerical  approximation  (equation  5 
below)  was  taken  from  Hanks  et  al . 
(1971). 

(Ti,j-Ti)j.l)/At=[(Ti.1J-TiJ)ai.1/2)j 
-(Ti,j-Ti+l,j)ai+l/2,j]/Az2-       [5] 


The  procedure  for  solving  equation  5  was 
similar  to  that  used  by  Hanks  et  al . 
(1971).   In  our  calculations,  AZ  was 
varied  from  small  values  near  the  soil 
surface  to  large  values  near  the  bottom 
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boundary.   Equation  5,  when  written  for 
each  depth,  results  in  a  tridiagonal 
matrix  that  was  solved  hourly  for  soil 
temperature  prediction  at  any  depth  and 
time . 

The  model  for  estimating  upper  boundary 
temperature  (UBT)  is  given  in  equation  6 


To,t=Ta  t+ATsin(27rtx/p>- 


6] 


Equation  6  is  the  summation  of  the  air 
temperature  at  a  reference  height  and  a 
sine  function  that  describe  the 
difference  in  the  temperature  of  air  and 
soil  surface  (fig.  6.1).   However,  this 
does  not  mean  that  the  soil  surface 
temperature  is  necessarily  a  sine 
function.   Presence  of  the  T 


a,  t 


term  in 


equation  6  implies  that  hourly 
temperature  fluctuations  at  the  soil 
surface  due  to  cloud  cover  or  wind  are 
reflected  in  the  air  temperature  at  the 
reference  height.   The  reference  height 
in  our  experiment  was  2  m.   Amplitude  of 
the  sine  function  (A<p)  in  equation  6  and 
figure  6.1  is  equal  to  the  difference 
between  maximum  temperatures  at  the  soil 
surface,  Tsx,  and  air,  Tax,  or  between 
the  minimum  temperatures  at  the  soil 
surface,  Tsn,  and  air,  Tan,  for  curves 
above  and  below  the  reference 
temperature,  respectively.   Period,  P,  is 
twice  the  hourly  difference  between  the 
time  that  the  temperatures  of  air  and 
soil  surface  are  similar.   A  reference 


SOIL  SURFACE 


21C° 


Figure  6.1. 

General   relationship  between  hourly  air 

and  soil   surface   temperature. 


temperature   of   21°C  was   assumed.       In 
equation   6 ,    P  would  have   different  values 
for  predicting   the   upper  boundary 
temperatures   above   and  below   the 
reference   temperature    for   each   day.      Thus 
in  figure   6.1, 


To,t=Ta,t+[Tsx-Taxsin(2*tx/P)] , 
P=2(t2-t1), 
and 


tx=t-tl 


[7] 
[8] 

[9] 


for  air  temperatures  above  the  reference 
temperature ,  and 

To,t=Ta,t+(Tsn-Tan)sin(27rtx/P),     [10] 

P=2(t3-t2),  and  [11] 

tx=t-t2  [12] 

for  air  temperature  below  the  reference 
temperature . 

Hourly  air  temperatures  were  estimated 
(Gupta  et  al.  1982)  from  the  daily 
maximum  and  minimum  air  temperatures 
using  equation  13. 


Ta,t=T?(Tax-Tan)+T 


an' 


13 


Hourly  "ft  values  are  taken  from  table 
6.1.   The  procedure  used  in  calculating 
t  is  given  in  the  next  section. 


T? 


Lower  boundary  conditions  are 
approximated  by  assuming  a  deep  profile 
such  that  temperatures  at  the  lower 
boundary  are  nearly  constant.   In  our 
simulation,  we  assumed  a  constant 
temperature  of  10°C  at  a  depth  of  5.5  m 
throughout  the  growing  season  and  for  all 
treatments.   This  assumption  was 
consistent  with  a  few  temperatures 
measured  over  the  growing  season  at  the 
2-m  depth  from  Morris,  MN. 
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Table  6.1. 

Averaged  normalized  hourly  air 
temperatures  (\t)    and  the  standard 
deviation  of  the  averaged  normalized 
hourly  air  temperatures  ("["£  D)  for  the 
growing  season  1978  and  197§  at  two 
locations  in  Minnesota  (J^±J^  n) 


Bulk  densities  ranged  from  1.2  g  cm"-5 


near  the  surface  to  1.59  g  cm' 
150 -cm  depth. 


at  the 


Morris 

.  MN 

St.  Paul,  MN 

Time 

1978 

1979 

1980 

0000 

0.25±0.19 

0.2610.21 

0.3610.28 

0100 

.21±.19 

.221.19 

.281.25 

0200 

.17±.17 

.191.17 

.231.23 

0300 

.14±.15 

.141.15 

.181.22 

0400 

.111.14 

.121.14 

.141.21 

0500 

.101.13 

.101.14 

.121.19 

0600 

.131.12 

.131.13 

.181.17 

0700 

.251.14 

.251.13 

.301.16 

0800 

.391.16 

.401.15 

.401.18 

0900 

.521.18 

.551.15 

.501.20 

1000 

.631.18 

.671.14 

.601.21 

1100 

.721.18 

. 741 . 14 

.681.21 

1200 

.761.16 

.831.13 

.781.18 

1300 

.831.14 

.881.13 

.851.17 

1400 

.881.14 

.921.10 

.881.19 

1500 

.911.14 

.931.11 

.911.18 

1600 

.921.14 

.921.13 

.891.18 

1700 

.881.15 

.881.15 

.841.20 

1800 

.791.18 

.811.16 

.771.20 

1900 

.661.19 

.651.19 

.661.20 

2000 

.491.17 

.501.19 

.531.21 

2100 

.381.19 

.401.20 

.441.23 

2200 

.321.20 

.331.21 

.371.25 

2300 

.261.20 

.271.23 

.311.26 

Verification 

Upper  boundary  and  root  zone  temperature 
predictions  were  tested  against  the 
measured  temperatures  from  the  1978 
residue  management  experiment  at  the  U.S. 
Department  of  Agriculture,  North  Central 
Soil  Conservation  Research  Center  at 
Morris,  MN .   The  soil  is  Hammerly  silty 
clay  loam  (Aerie  Calciaquolls) .   The 
surface  soil  is  18.4  percent  sand,  53.0 
percent  silt,  28.8  percent  clay,  and  6.2 
percent  organic  matter.   Subsurface 
samples  varied  widely  in  contents  of  sand 
(20  to  81  percent),  silt  (14  to  46 
percent),  and  clay  (14  to  38  percent). 


The  experimental  site  was  tilled  May  11, 
1978,  and  corn  planted  on  May  12,  1978. 
The  oat  (Avena  sativa)  straw  mulch  plots 
were  established  on  May  17,  1978. 

The  residue  experiment  consisted  of  six 
treatments  representing  various  degrees 
of  cover  resulting  from  either  residue  or 
shading  by  growing  plants.   The 
treatments  were  a  bare  soil  surface 
without  a  crop,  a  50 -percent  oat  straw 
mulch  cover  without  a  crop,  a  75 -percent 
oat  straw  mulch  cover  without  a  crop,  a 
99 -percent  oat  straw  mulch  cover  without 
a  crop,  a  bare  soil  surface  with  a  corn 
(Zea  mays  L. )  crop,  and  a  50 -percent  oat 
straw  mulch  cover  with  a  corn  crop. 

The  amount  of  straw  was  estimated  from 
the  plot  area  and  the  relationship  given 
by  Sloneker  and  Moldenhauer  (1977) : 

percent  cover=100(l-e   '     ),      [14] 

where  x  is  the  rate  of  oat  straw  mulch  in 
metric  tons  per  hectare.   The  plots  were 
15  m  long  and  9  m  wide.   Row  spacing  for 
corn  treatments  was  76  cm  with  the  rows 
running  east  to  west.   Corn  variety  was 
Pioneer  3955,  and  the  population  was 
56,790  plants  per  hectare. 

Soil,  water,  and  temperatures  were 
monitored  throughout  the  growing  season. 
Copper- cons tantan  thermocouples  were 
buried  at  the  surface  and  5,  10,  15,  20, 
and  30  cm  deep  in  all  plots .   Four 
parallel  thermocouples,  randomly  located 
in  the  plot,  spatially  averaged  the  soil 
temperature  at  each  depth.   Neutron  probe 
measurements  (down  to  150  cm)  were  taken 
once  a  week,  and  moisture  block  and 
tensiometer  readings  were  taken  three 
times  a  week.   Gravimetric  soil  water 
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content  was  also  determined  several  times 
during  the  season.   Particle  size 
analysis,  organic  matter  content,  and 
bulk  densities  were  used  to  estimate 
soil  water  retention  characteristic 
curves  at  each  depth  (Gupta  and  Larson 
1979). 

Thermal  dif fusivities  were  calculated 
from  the  soil  mineralogy,  organic  matter 
content,  and  bulk  density  using  DeVries' 
(1963)  relationship  and  modifications 
(curves  A  and  B  in  figure  1  of  Kimball  et 
al . ,  1976).   In  the  model,  the  thermal 
diffusivities  varied  with  depth  but  were 
kept  constant  over  the  growing  season. 
Thermal  diffusivities  at  each  depth  cor- 
respond to  an  average  water  content  over 
the  growing  season  as  estimated  from  the 
tensiometer  readings  and  water  retention 
curves.   Future  versions  of  the  soil 
temperature  submodel  could  be  run  in 
conjunction  with  the  NTRM  main  model  and 
could  adjust  a  as  a  function  of  water 
content. 

We  used  an  average  matric  potential  of 
-0.2  bar  at  15 -cm  depth  and  0  at  150 -cm 
depth  for  the  calculation  of  thermal 
diffusivities.   However,  there  were 
several  days  during  the  growing  season 
when  the  matric  potential  was  outside  the 
measurement  range  of  tensiometers  (<-0.8 
bar)  in  the  top  70  cm  of  the  profile.   In 
the  top  2  cm,  we  used  thermal  diffusivi- 
ties corresponding  to  an  air-dry  water 
content  (0=0.07).   In  addition,  it  was 
assumed  that  water  content  between  depths 
of  2  and  15  cm  corresponded  to  -0.2  bar 
matric  potential  and  that  below  the  150- 
cm  depth,  the  soil  profile  was  homoge- 
neous, with  a  thermal  diffusivity 
corresponding  to  that  at  the  150 -cm 
depth.   The  same  thermal  diffusivity 
distribution  with  soil  depth  was  used  for 
all  treatments  throughout  the  growing 
season. 

All  micrometeorological  data,  including 
air  temperature  at  several  heights,  were 
gathered  at  a  weather  station  immediately 


next  to  the  plots.   The  soil  surface  at 
this  weather  station  was  grass  covered. 
Hourly  air  temperatures  (2-m  height) 
measured  at  this  weather  station  were 
used  in  the  calculation  of  soil  upper 
boundary  temperatures  for  all  treatments . 
This  assumes  that  air  at  2-m  height  was 
well  mixed  over  the  plot  area  and  was  not 
affected  by  residue  or  crop  cover  at  the 
soil  surface.   All  the  micrometeoro- 
logical data  and  soil  temperatures  were 
logged  either  hourly  or  semihourly  by 
automatic  data  acquisition  systems  from 
June  13  through  September  24. 

Hourly  air  temperatures  were  normalized 
with  respect  to  daily  maximum  and  minimum 
air  temperature  using  equation  15. 


Ta, t~(Ta, t"Tan)/(Tax"Tan) ' 


15 


where  Ta  t  =  t^ie  normalized  hourly  air 

temperature, 
Ta  t  =  the  hourly  air  temperature, 

°C, 
Tan  =  the  daily  minimum  air 

temperature,  °C,  and 
Tax  =  the  daily  maximum  air 

temperature,  °C. 

Normalized  hourly  air  temperatures  were 
then  averaged  over  the  growing  season 
using  equation  16. 

T?=ITa,t/N  [16 

i=l 

where  ]"£  =  the  normalized  air 

temperature  at  t  hour 
averaged  over  the  growing 
season  and 
N  -  the  number  of  observations . 

Standard  deviation  of  the  averaged 
normalized  hourly  air  temperature  (Tt,D^ 
was  calculated  from  equation  17 . 


T?,D=J  [l(Ta,t)2-NT?2]/(N-D.      .  [17] 
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Figure  6.2  is  the  plot  of  ]~*  at  each  hour 
of  the  day  for  three  growing  seasons  at 
two  locations.   The  data  in  figure  6.2 
cover  a  period  of  May  25  through 
September  27,  1978  (curve  A),  and  June  13 
through  September  24,  1979  (curve  B) ,  at 
Morris,  MN,  and  June  24  through  September 
22,  1980  (curve  C) ,  at  St.  Paul,  MN. 
Both  table  6.1  and  figure  6.2  show  that 
t  values  were  about  the  same  for  these 
growing  seasons  and  locations.   Minor 
differences  in  the  "[2  values  occurred 
between  Morris  and  St.  Paul  in  the  early 
and  late  hours  of  the  day.   Figure  6.2 
shows  that  the  averaged  maximum 
temperatures  occur  at  1600  hours  for  the 
growing  season  of  1978  at  Morris  and  1500 
hours  for  the  growing  season  of  1979  and 
1980  at  Morris  and  St.  Paul, 
respectively.   Averaged  minimum 
temperatures  occurred  at  0500  hours  for 
all  three  growing  seasons  from  the  two 
locations.   At  any  hour  of  the  day,  Jt   n 
were  generally  lower  for  Morris  than  St. 
Paul  (table  6.1).   Also  for  a  growing 
season,  Jt  p  decreased  and  increased 
during  the  day  with  the  minimum  values 
occurring  near  the  time  of  daily  maximum 
air  temperature. 


6/13/78  -  9/2a/78  Moms.  MN 

6/24/79  -  9/22/79  Moms.  MN 

5/25/80  -  9/27/60  St   Paul,  MN 


0800 


TIME,  HRS. 
Figure   6.2. 

Averaged  normalized  hourly  air 
temperatures  (j£  for  the  growing  seasons 
of  1978  and  1979  at  Morris,  MN,  and  1980 
at  St.  Paul,  MN) . 


Similarity  between  the  three  curves 
(figure  6.2)  from  Morris  and  St.  Paul 
suggests  that  even  for  a  broader  region 
like  the  northern  Corn  Belt,  the  hourly 
air  temperatures  when  normalized  between 
the  daily  maximums  and  minimums  may 
result  in  curves  similar  to  curves  in 
figure  6.2. 

Periods,  P,  were  estimated  by  selecting 
*-l '  t2 '  anc*  t3  (figure  6.1)  such  that  the 
estimated  hourly  air  temperatures  were 
either  equal  to  or  slightly  higher  than 
the  reference  temperature . 


Residue  Cover  vs.  Maximum  and  Minimum 
Surface  Temperature 

Figure  6.3  shows  the  relationship  between 
the  daily  maximum  and  minimum 
temperatures  of  air  (at  2-m  height)  and 
of  residue -covered  soil  surfaces.   Solid 
lines  were  obtained  by  fitting  a  cubic 
spline  function  to  the  measured  values. 
From  June  13  through  September  24, 
maximum  air  temperatures  ranged  from 
37.5CC  to  14.4CC,  whereas  minimum  air 
temperatures  varied  between  1.0°C  to 
21.4°C.   Figure  6.3  is  shown  in  a  split 
graph  only  to  amplify  the  differences 
between  the  treatments .   For  a  given  air 
temperature,  the  maximum  temperature  of 
the  soil  surface  (figure  6.3b)  decreased 
and  the  minimum  temperature  (figure  6.3a) 
of  the  soil  surface  increased  with 
increasing  amounts  of  soil  surface  cover. 
This  is  similar  to  the  observations  of 
Van  Wijk  et  al.  (1959)  and  Burrows  and 
Larson  (1962) . 

Lower  maximum  temperatures  are  the  result 
of  lower  thermal  diffusivity  of  surface 
residues  and  greater  reflection  of  solar 
radiation  with  increase  in  plant  or 
residue  cover.   Higher  minimum 
temperatures  are  the  result  of  a  decrease 
in  the  long  wave  radiation  emission  at 
night  from  plant-or  residue-covered  soil 
surfaces.   Average  heights  of  a  corn 


30 


plant  were  22  cm,  218  cm,  and  241  cm  on 
May  30,  July  15, and  August  6, 
respectively.   The  switch  between  curves 
2  and  3  at  air  temperatures  greater  than 
28°C  is  a  result  of  large  scatter  in  the 
measured  temperature  for  the  75 -percent 
residue  covered  treatment  (curve  3) .   The 
maximum  (Tgx)  and  the  minimum  (Tsn) 
temperature  of  the  soil  surface  for  a 
desired  residue-cover  treatment  were 


estimated  from  the  maximum  (Tax)  and 
minimum  (Tan)  air  temperatures  and  the 
curve  in  figure  6.3. 


Hourly  Air  Temperature  Predictions 

Figure  6.4  shows  the  comparison  between 
the  predicted  and  measured  hourly  air 
temperatures  at  2-m  height  over  a  sod 


AIR  TEMPERATURE  (°C)  AT  2m 


Figure  6.3. 

Relationship  between  the  daily  maximum  (a) 
and  minimum  (b)  temperatures  of  air  (20-m 
height)  and  of  bare  and  residue-covered 
soil  surfaces  with  and  without  a  growing 
corn  crop. 
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surface,  for  a  4-day  period.   Predicted 
hourly  air  temperatures  are  discontinuous 
from  one  day  to  the  next.   This 
discontinuity  arises  because  the  maximum 
and  the  minimum  air  temperatures  are 
different  in  two  consecutive  days. 
Maximum  and  minimum  predicted  air 
temperatures  (fig.  6.4)  are  lower  and 
higher  than  the  measured  maximum  and 
minimum  air  temperatures  (input), 
respectively.   These  differences  result 
because  the  values  of  ~|~£  (curve  A  in 
figure  6.2  and  table  6.1)  corresponding 
to  minimum  and  maximum  air  temperature 
are  0.1  and  0.92,  respectively.   Ideally, 
these  values  would  be  0.0  and  1.0,  if  the 
minimum  and  maximum  air  temperatures 
occurred  at  the  same  time  every  day 
during  a  season.   Differences  between  the 
measured  and  predicted  air  temperatures 
in  figure  6.4  are  less  than  4°C. 


Since  air  temperatures  are  also  used  to 
calculate  the  growing  degree  days  in  crop 
growth  modeling,  comparisons  were  also 
made  between  the  values  for  growing 
degree  days  calculated  from  both  the 
measured  and  predicted  air  temperatures 
(table  6.2).   Growing  degree  days  were 
calculated  using  equations  20  and  21. 


GDD=X(Ta  t-10.0)/24.0 


for  10<Ta  t<35°C. 


20] 


GDD=X[((Tax+Tan)/2.0)-10.0].        [21] 

Table  6.2  shows  that  the  errors  in  the 
calculation  of  growing  degree  days  are 
less  than  5  percent  when  measured  air 
temperatures  were  replaced  with  the 
predicted  air  temperature  in  equations  20 
and  21,  respectively. 


Closeness  of  the  predicted  temperatures 
to  the  measured  values  is  also  checked  by 
comparing  the  difference  between  the 
predicted  and  measured  temperatures 
averaged  over  the  growing  season  with  the 
standard  deviation  of  the  average 
difference.   The  average  difference  (D) 
is  defined  as, 


D=l/N£(Tp-Tm), 


18 


and  the  standard  deviation  of  the  average 
difference  (Sp)  is  defined  as 


■>^n 


(Tp-Tm) 


Li' 


•ND 


-]/(N-l), 
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where  Tp  and  Tm  are  the  predicted  and 
measured  temperatures,  respectively,  and 
N  is  the  number  of  observations . 

The  average  difference  and  the  standard 
deviation  of  the  average  difference  for 
the  hourly  air  temperature  are  0.08°C  and 
1.96CC,  respectively,  over  the  growing 
season.   Thus,  95  percent  (D±1.96SD)  of 
the  predicted  hourly  air  temperatures  are 
within  -3.8°C  to  3.9°C  of  the  measured 
values . 


Daily  Root  Zone  Temperature  Predictions 

Figures  6 . 5  through  6 . 7  show  the 
difference  between  the  daily  predicted 
and  measured  root  zone  temperatures  for 
bare,  50-percent-oat-residue-cover  and 
corn-plus -50 -percent -oat -residue -cover 
soils.   Root  zone  temperatures  are 
predicted  from  using  curve  A  of  figure 
6.2.   During  the  growing  season,  root 
zone  temperatures  predicted  from  the 
temperature  model  followed  the  general 
trend  of  the  measured  values . 

Table  6 . 3  gives  the  value  of  D  and  SD  of 
the  root  zone  temperatures  when  measured 
and  estimated  upper  boundary  temperatures 
are  used  in  the  numerical  analysis  of  the 
diffusion  equation  for  various  residue 
cover  treatments .   Data  for  measured 
upper  boundary  temperatures  in  table  6.3 
are  taken  from  Gupta  et  al .  (1981).   When 
measured  upper  boundary  temperatures  are 
used  in  the  numerical  solution  of 
diffusion  equation,  95  percent  (D±1.96 
SD)  of  the  predicted  daily  temperatures 
at  5-cm  depth  are  within  -1.3°C  to  3.0°C 
and  -1.1°C  to  1.2°C  of  the  measured 
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Table  6.2. 

Growing  degree  days  as  calculated  from  the 

measured  and  predicted  air  temperatures 


Method  of  calculating 
growing  degree  days 


170 


190 


Julian  day 


210 


230 


250 


267 


Equation  20: 
Measured 
Predicted 

Equation  21: 
Measured 
Predicted 


70 

290 

513 

729 

957 

1,062 

66 

285 

501 

712 

953 

1,060 

68 

289 

508 

722 

973 

1,075 

68 

291 

513 

730 

983 

1,089 

Table  6.3. 

Average  difference  (D)  and  standard  deviation  of  the 
average  difference  (Sn)  between  daily  predicted  and 
measured  root  zone  temperatures  at  5- ,  10-  and  30 -cm 
depths  as  influenced  by  measured  or  estimated 
(simplified  method)  upper  boundary  temperatures  for  bare 
and  residue-covered  soil  surfaces  with  and  without  a 
growing  corn  crop 


Methods  of 

determining 

upper 

boundary 

temperature 


D±SD  (°C) 


Bare  soil   50%  Cover   75%  Cover   99%  Cover 


Corn,  no 
residue 


Corn,  50^ 
residue 


5 -cm  depth 


Measured 
Estimated 


0.85±1.10 
.90±3.16 


0.47±0.81 
. 20±2 . 54 


1.01±1.27 
1.43±2.92 


0.09±0.32 
.4911.36 


•0.1010.50 
.2811.29 


0.0610.57 
.061  .12 


10-cm  depth 


Measured 
Estimated 


0.7110.88 
.7612.60 


0.6210.59 
.3512.07 


0.8310.85 
1.2212.26 


0.2010.24 
.5711.08 


0.0810.36 
.4211.22 


0.0010.42 
.3511.20 


30 -cm  depth 


Measured 
Estimated 


1.0410.51 
1.0711.57 


0.6610.37 
.3811.30 


1.2110.52 
1.5211.35 


0.6310.30 
.911  .69 


0.8510.39 
1.081  .98 


0.1510.23 
.091  .92 
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TIME,  HRS. 

Figure  6.4. 

Comparison  between  the  measured  and 
predicted  hourly  air  temperatures  for  a 
4-day  period. 

values  for  bare  and  corn-with-50-percent- 
residue-cover  soils  respectively.   These 
ranges  increase  by  4°C  for  bare  soil  and 
by  2°C  for  corn-with-50-percent-residue- 
cover  soil  when  measured  upper  boundary 
temperatures  are  replaced  with  boundary 
temperatures  estimated  from  the  procedure 
described  in  this  report  model.   Thus,  95 
percent  of  the  daily  temperatures  at  5-cm 
depth  predicted  using  the  upper  boundary 
temperature  estimated  from  the 
simplified  model  are  within  -5.3°C  to 
7.1°C  and  -2.0°C  to  3.0°C  of  the  measured 
value  for  bare  and  corn-with-50-percent- 
residue-cover  soils,  respectively. 

In  general,  the  difference  between 
predicted  and  measured  root  zone 
temperatures  (table  6.3)  decreases  with 
an  increase  in  (1)  the  amount  of  residue 
cover  at  the  soil  surface  and  (2)  the 
soil  depth. 

Sensitivity  Analysis 

Table  6.4  gives  the  sensitivity  analysis 
of  the  root  zone  temperature  predictions 
.nfluenced  by  the  difference  in  the 
.  alized  air  temperature  curves  (fig. 
6.2)  from  two  locations  and  three  growing 
seasons .   The  results  are  presented  in 
terms  of  D  and  Sjj  for  the  bare  soil. 


204  214  224 

JULIAN  DAY 

Figure   6.5. 

Daily  measured  and  predicted  soil 
temperatures  at  three  depths  for  bare 
soil  during  the  growing  season.   Upper 
boundary  temperatures  used  in  the 
numerical  analysis  were  estimated  by  the 
simplified  method. 


o-     10 
—     30 


I  2°C 


214  224 

JULIAN  DAY 


Figure   6.6. 

Daily  measured  and  predicted  soil 
temperatures  at  three  depths  for  50- 
percent  residue -covered  surface  during 
the  growing  season.   Upper  boundary 
temperatures  used  in  the  numerical 
analysis  were  estimated  by  the  simplified 
method.- 
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204  214         224 

JULIAN  DAY 
Figure  6.7. 

Daily  measured  and  predicted  soil 
temperatures  at  three  depths  for  50- 
percent  residue -covered  surface  plus  corn 
during  the  growing  season.   Upper 
boundary  temperatures  used  in  the 
numerical  analysis  were  estimated  by  the 
simplified  method. 

Table  6.4  shows  that  the  variations  in 
the  predicted  root  zone  temperatures  are 
slight  when  upper  boundary  temperatures 
estimated  from  curve  A,  curve  B,  or  curve 
C  (fig.  6.2)  are  used  in  the  numerical 
analysis  of  the  diffusion  equation.   This 
indicates  that  any  of  the  normalized  air 
temperature  curves  in  figure  6.2  may  be 

Table  6.4. 

Average  difference  (D)  and  standard 
deviation  of  the  average  difference  (Sp) 
between  the  daily  predicted  and  measured 
root  zone  temperatures  of  a  bare  soil  at 
5-,  10-  and  30 -cm  depth  as  influenced  by 
the  estimates  of  upper  boundary 
temperatures  calculated  using  curves  A, 
B,  and  C  of  figure  6.1 


used  to  predict  the  root  zone  temperature 
at  other  locations  without  causing 
significant  additional  errors.   The  input 
variables  that  would  vary  with  location 
are  daily  maximum  and  minimum  air 
temperatures  and  thermal  dif fusivities  of 
the  soil  profile. 


Summary  and  Conclusions 

We  describe  models  that  can  be  used  to 
predict  root  zone  temperature  of  bare  and 
residue -covered  soil  surfaces  with  and 
without  a  corn  crop.   Input  to  the  model 
are  daily  maximum  and  minimum  air 
temperatures  at  a  2-m  height,  initial 
soil  temperature  profile,  and  thermal 
diffusivity  variation  with  depth. 
Thermal  dif fusivities  are  estimated  from 
soil  mineralogy,  organic  matter  content, 
bulk  density,  and  field  measurements  of 
matric  potential. 

Root  zone  temperatures  predicted  from  the 
soil  temperature  model  followed  the 
general  trend  of  the  measured  values 
during  the  growing  season.   Ninety- five 
percent  of  the  daily  predicted 
temperatures  at  5 -cm  depth  are  within 
-5.3°C  to  7.3°C  and  -2.0°C  to  3.0°C  of 


Depth 
(cm) 

Curve  A 
(D+SD  (°C)) 

Curve  B 
(D+SD  CO) 

Curve  C 
(D+SD  CO) 

5 
10 
30 

0.90±3.16 

.76±2.60 

1.07±1.57 

1.08±3.15 

.93±2.58 

1.22±1.56 

1.2213.30 

1.0712.70 
1.3411.63 
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the  measured  values  for  bare  and  corn- 
plus-50-percent-residue-cover  soils , 
respectively.   Differences  between  the 
predicted  and  measured  root  zone 
temperature  decreased  with  an  increase  in 
(1)  plant  or  residue  cover  at  the  surface 
and  (2)  soil  depth. 

The  soil  temperature  model  can  be  a 
valuable  tool  in  predicting  spring  root 
zone  temperatures  for  various  plant  or 
residue  covers  in  the  Northern  Corn  Belt 
States.   Seven- day  forecasts  of  the  daily 
maximum  and  minimum  air  temperatures  can 
be  used  in  the  soil  temperature  model  to 
provide  an  approximation  of  the  daily 
root  zone  temperature  under  various 
residue  covers.   Daily  root  zone 
temperature  prediction  can  then  be 
continuously  improved  over  the  growing 
season  by  replacing  the  forecast  for  the 
previous  week  with  the  measured  values  of 
daily  maximum  and  minimum  air 
temperature.   Forecasting  of  the  daily 
root  zone  temperatures  under  various 
residue  covers  should  be  useful  in 
deciding  the  planting  dates  for  the 
Northern  Corn  Belt  States. 
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Chapter  7 

CARBON  AND  NITROGEN  TRANSFORMATIONS  IN 

SOIL,  SUBMODEL  NCSOIL 

J.A.E.  Molina,  C.E.  Clapp,  M.J.  Shaffer, 
F.W.  Chichester,  and  W.E.  Larson1 

NCSOIL  is  the  submodel  of  NTRM  that 
describes  the  carbon  and  nitrogen 
transformations  in  soil.   NCSOIL  computes 
the  changes  in  carbon  and  nitrogen 
organics,  ammonium  (exchangeable  plus 
soluble) ,  and  nitrate  concentrations  that 
result  from  the  processes  of  residue 
decomposition,  mineralization, 
immobilization,  nitrification, 
denitrif ication,  and  nonsymbiotic 
nitrogen  fixation.   The  rates  depend  on 
soil  water  and  temperature.   A  mechanism 
is  included  to  account  for  the  increase 
in  the  rate  of  mineralization  after  a 
physical  disturbance  (tillage)  of  the 
soil.   NCSOIL,  which  was  devised  for 
short-term  N  and  C  dynamics,  is  a 
theoretical  model  (Weigert  1975)  built  on 
the  concept  of  catenary  sequence  of 
heterogeneous  substrates  (Van  Veen  et  al . 
1981) .   The  soil  organic  substrates  are 
defined  by  their  kinetic  constants  and 
their  position  in  the  model's  structure. 
Residues  are  defined  in  reference  to 
entities  that  are  chemically  or 
morphologically  static. 

Input  variables  and  parameters  internal 
to  NCSOIL  were  kept  to  a  minimum, 
compatible  with  the  necessity  to  maintain 
within  NCSOIL  a  mechanistic  structure  to 


compute  a  variable  efficiency  factor  and 
total  as  well  as  isotopic  nitrogen 
kinetic  rates .   Seven  carbon  and  nitrogen 
pools  were  considered:   the  organic 
residues,  pool  I  and  pool  II  of  the  soil 
active  organic  phase,  exchangeable  plus 
soluble  ammonium,  nitrate,  dinitrogen, 
and  carbon  dioxide.   The  passive  pool  of 
soil  organics  (stable  humus)  was  assumed 
to  have  no  influence  on  the  dynamics  of 
the  system  over  the  short  period,  a 
growing  season,  considered  for  each  run 
of  the  model.   Some  components  of  NCSOIL 
are  similar  to  those  introduced  by  Beek 
and  Frisel  (1973)  and  Van  Veen  (1977). 
When  appropriate,  the  same  variable 
symbols  have  been  used.   NCSOIL  is 
written  in  FORTRAN  77,  and  summation  of 
the  rates  is  performed  by  a  Runge-Gutta' s 
method. 

Chichester  et  al .  (1975)  reported  the 
results  of  an  experiment  performed  to 
investigate  the  "relative  mineralization 
rates  of  indigenous  and  recently 
incorporated  15-N  labeled  nitrogen." 
Results  reported  in  Chichester  et  al . 
(1975)  as  well  as  unpublished  data  of  the 
same  authors  (in  particular,  rates  of  CO2 
released  during  the  preincubation 
periods)  were  used  to  analyze  the 
behavior  of  NCSOIL. 


1  Molina  is  a  professor,  Soil 
Science  Department,  University 
of  Minnesota,  St.  Paul,  MN 
55108;  Clapp  is  a  research 
chemist  and  Shaffer  is  a 
research  soil  scientist,  U.S. 
Department  of  Agriculture, 
Agricultural  Research  Service, 
Soil  and  Water  Management 
Research  Unit,  St.  Paul,  MN 
55108;  Chichester  is  a 
research  soil  scientist,  U.S. 
Department  of  Agriculture, 
Agricultural  Research  Service, 
Temple,  TX  76503;  and  Larson 
is  head,  Soil  Science 
Department,  University  of 
Minnesota,  St.  Paul,  MN  55108. 


Model  Description 

The  active  organic  phase  (Jansson  1958) 
was  divided  into  two  substrates,  pool  I 
and  pool  II.   Pool  I  serves  as  a  sink  for 
the  organic  decomposition  products  of  the 
residues  and  for  a  fraction  of  the 
decomposition  products  of  pool  I  and  pool 
II  (figs.  7.1-7.3).   Pool  II  is  the 
recipient  of  that  organic  fraction  of 
pool  I  that  is  not  recycled.   A  single 
pool  I  is  considered:   no  distinction  is 
made  among  the  decomposition  constants  of 
the  various  pools  that  may  arise  from  the 
decomposition  of  various  substrates.   It 
was  assumed  that  the  passive  (resistant) 
organic  phase  (stable  humus)  would  not 
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interact  with  the  dynamics  of  the  system 
within  the  limited  time,  a  growing 
season,  during  which  the  model  is  run. 
Therefore  EFFHUM  =0.0,  and  only  three 
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Figure  7.1. 

Flowchart  for  nitrogen  transformations. 
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Figure   7.2. 

Relational   diagram   for   the   decomposition 
of  a  carbon  residue   component    (labile   or 
resistant) .      The   temperature   and  water 
reduction  factors   are  not   shown.      The 
symbols   correspond  to   those  used   in 
system-dynamic   simulation   (Beek  and 
Frisel   1973,    Ferrari  1978)    and  are 
explained  briefly   in  the   text. 


organic  substrates  are  available  for 
microbial  degradation:  the  residues, 
pool   I,    and  pool    II. 

Each  heterogeneous    substrate  was   handled 
as    if   two   components,    one   labile   and   the 
other  resistant,    were   decomposing 
exponentially  and   independently   from  each 
other    (Eckenfelder   1970,    Hunt   1977,    and 
Jenkinson   1977).      NCSOIL   in  NTRM  has 
provisions   for   20  residues    (40 
components) ;    as   a   stand-alone   program,    it 
provides   for  only  one   residue.      Thus   the 
potential   rate   of  decomposition  of  a 
substrate   is   given  by: 

dC/dt(potential)=/i[S»k»C+(l-S)h.C]  ,  [1] 


where 


C  = 


S  and  (1-S) 


the  carbon  concentration  of 
the  substrate  at  the  time 
t, 

the  time  invariant 
proportions  of  the  labile 
and  resistant  components  of 
the  substrate, 
k  and  h  =  the  specific  rates  of 
decomposition  in  each 
component , 
and  /i  =  the  lowest  of  the 

temperature  or  water 
reduction  factors. 

The  difference  between  the  actual  and 
potential  rates  of  residue  decomposition 
(fig.  7.1)  reflects  a  reduction  in  the 
rate  of  decomposition  when  not  enough 
nitrogen  is  available.   Therefore, 

dC/dt( actual )=^NdC/dt (potential) ,    [2] 

where  /j,-^   is  a  reduction  factor  that  is  a 
function  of  the  ratio  of  the  potential 
rate  of  carbon  residue  decomposition  to 
the  total  amount  of  nitrogen  available  (N 
potentially  released  from  the  residue 
during  the  decomposition  plus  the 
inorganic  nitrogen  available) .   The 
variable  PX  (fig.  7.2)  is  identical  to 
the  one  introduced  by  Beek  and  Frisel 
(1973)  to  determine  whether  residue 
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Relational  diagram  for  carbon  pool  I  and 

pool  II  decomposition.   The  temperature 

and  water  reduction  factors  are  not 

shown. 
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decomposition  will  result  in  nitrogen 
mineralization  or  immobilization.   In  a 
situation  of  nitrogen  immobilization 
(PX<0) ,  the  surplus  nitrogen  that  has  to 
be  obtained  from  the  inorganic  pool  or 
atmospheric  nitrogen  is  computed  as 
SURPN.   This  amount  is  a  function  of  the 
maximum  efficiency  factor  (EFFAC) ,  which 
is  reduced  to  an  actual  efficiency  factor 
to  balance  SURPN  with  the  inorganic 
nitrogen  actually  available  or  the  rate 
of  nonsymbiotic  nitrogen  fixation  when 
the  inorganic  pool  is  depleted.   In 
NCSOIL,  the  efficiency  factor  associated 
with  residue  decomposition  is  variable. 

At  each  computational  time  step,  the 
carbon  recycled  into  pool  I  is  divided 
into  its  labile  and  resistant  portions. 
In  contrast,  the  fraction  (EFFNO)  of 
decomposed  pool  I  feeds  directly  into  the 
resistant  pool  II.   Every  time  the  soil 
is  physically  disturbed  (tillage  or  soil 
sampling  and  preparation  for  laboratory 
incubation) ,  the  tillage  input  flag 
(TILL)  is  raised  and  a  time  invariant 
fraction  of  the  resistant  pool  II  is 
transferred  to  the  labile  pool  II;  thus, 
the  rate  of  mineralization  is  increased. 
Since  the  C:N  ratio  of  pool  II  is  chosen 
equal  to  the  C:N  ratio  of  pool  I,  the 
decomposition  of  pool  I  and  of  pool  II 
results  in  nitrogen  mineralization  (Beek 
and  Frisel  1973) ,  and  their  corresponding 
efficiency  factors  are  set  to  the  maximum 
value  EFFBAC.   It  also  follows  that 
immobilized  nitrogen  can  only  merge  with 
the  flow  of  organics  between  the  residues 
and  pool  I  (fig.  7.3).   The  total 
inorganic  nitrogen  demand  due  to 
immobilization  is  distributed  over  the 
ammonium  and  nitrate  pools  in  proportion 
to  their  respective  concentrations. 

For  a  given  soil,  the  rate  of 
nitrification  is  constant  (Sabey  et  al . 
1969) .   When  the  water  content  of  the 
soil  has  reached  saturation, 
nitrification  is  replaced  by  the  process 
of  denitrif ication  with  a  rate 
proportional  to  the  amount  of  energy 


available,  as  measured  by  the  sum  of  the 
rates  of  decomposition  of  carbon 
residues,  pool  I,  and  pool  II. 
Nitrification  and  denitrif ication  can  be 
simultaneously  activated.   Boundary 
conditions  arise  when  immobilization 
competes  with  nitrification  for  ammonium 
or  with  denitrif ication  for  nitrate.   The 
model  assumes  dominance  of  immobilization 
over  nitrification  for  ammonium  uptake. 
When  the  initial  nitrate  supply  is  not 
sufficient  to  satisfy  both 
denitrif ication  and  immobilization,  the 
carbon  flow  rates  are  recomputed  with  a 
nitrate  level  reduced  by  the  nitrate 
demand  for  denitrif ication.   After  this 
correction,  precedence  is  given  to 
immobilization  over  denitrif ication. 
Nonsymbiotic  nitrogen  fixation  is 
activated  only  when  both  immobilization 
calls  for  ammonium  and  nitrate  and  the 
inorganic  nitrogen  pool  is  low.   The  rate 
of  fixation  is  proportional  to  the  amount 
of  energy  available,  as  measured  by  the 
potential  rate  of  carbon  residues 
decomposition. 

Rates  of  tagged  N  transformations  are 
computed  by  multiplying  the  rates  of  the 
corresponding  total  N  transformations  by 
the  ratio  of  the  tagged  N  to  total  N 
concentrations  of  the  substrates. 


Calibration  and  Behavior  of  NCSOIL 

Pools  I  and  II  are  given  operational  and 
dynamic  definitions  by  their  kinetic 
constants  and  their  position  in  the  C  and 
N  flowcharts.   The  active  organic  phase 
is  complex;  some  attempts  have  been  made 
to  categorize  it  broadly  in  descriptive 
terms  based  on  chemical  properties  or 
morphological  differences  (McGill  et  al. 
1974).   It  comprises  micro-organisms, 
microbial  metabolites,  lytic  products  of 
dead  cells,  and  the  humads  (Eckenfelder 
1970).   Obviously,  pool  I  and  pool  II  are 
an  operational  simplification  of  the  soil 
active  organic  phase.   Yet  analogies,  if 
not  absolute  correspondence,  between  the 
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dynamic  substrate  pools  of  NCSOIL  and  the 
static-chemical-morphological  entities  of 
the  active  organic  phase  are  important  to 
delineate  for  any  tentative  comparison  of 
the  model's  decomposition  constants  with 
those  found  in  the  literature  and 
accruing  to  chemicals  or  morphological 
entities.   Pool  I  and  pool  II,  in  view  of 
their  position  in  the  model's  structure, 
can  be  viewed  as  the  soil  biomass  and  the 
humads ,  terms  sufficiently  vague  to  suit 
the  analogy  and  convey  the  quality  of  the 
model . 

The  distinction  within  each  substrate  of 
a  labile  and  a  resistant  component 
results  from  the  choice  of  equation  1. 
Should  the  rate  equation  be  composed  of 
the  sum  of  three  exponential  terms,  three 
components  for  each  of  the  residues,  pool 
I,  and  pool  II  would  have  to  be 
considered.   The  concentrations  of  the 
labile  and  resistant  components  are 
defined  by  equation  1  and  do  not  have  any 
material  support.   Those  concentrations 
reflect  a  dynamic  property  of  the  active 
organic  phase  in  the  same  way  as  the 
dynamic  mass  of  an  object,  obtained 
indirectly  by  dynamic  measurements ,  is  a 
different  concept  from  its  static  mass, 
which  can  be  measured  by  direct 
comparison  with  others  of  the  same  kind. 
Rate  equation  1  expresses  a  zero-order 
kinetics  with  respect  to  the  active 
biomass  (which  could  be  expressed  by  a 
constant  or  variable  fraction  of  pool  I) . 
The  influence  of  the  degradative  agent 
concentration  is  thus  relegated  in  the 
decomposition  constants,  k,  h,  and  S. 

The  exclusion  from  the  rate  equation  of 
the  biomass  as  an  independent  variable 
allowed  an  integral  function  to  be 
obtained- -concentration  =  f  (k,  h, 
S/time) - -ready  to  be  fitted  against  the 
numerous  data  about  substrate 
degradation,  which  are  presently 
available  in  the  literature  and  which 
usually  do  not  include  any  information 
about  the  biomass. 


The  decomposition  in  soil  of  whole  cells 
obtained  from  axenic  microbial  cultures 
is  rapid  and  exhibits  a  remarkable 
constant  kinetics  whether  bacteria, 
yeast,  or  fungi  are  considered  (Nelson  et 
al.  1979).   As  a  basis  for  the  calibra- 
tion of  the  model,  the  decomposition 
constants  (k,  h,  and  S)  for  whole  cells 
were  assigned  to  pool  I,  even  though,  to 
repeat,  it  does  not  follow  that  pool  I 
equates  with  whole  microbial  cells  from  a 
chemical  and  morphological  standpoint. 
Stanford's  potentially  mineralizable 
nitrogen  NQ  (Stanford  and  Smith  1972)  is 
a  dynamic  concept,  and  the  nitrogen 
content  of  pool  I  plus  pool  II  does 
correspond  to  NQ.   The  constants  corre- 
sponding to  the  decomposition  of  NQ 
appear  to  be  independent  of  the  soil 
type,  provided  that  the  net  mineraliza- 
tion kinetics  were  obtained  from  soils 
that  had  not  been  recently  enriched  with 
an  exogenous  residue  input  (Stanford  and 
Smith  1972) ,  a  situation  corresponding  to 
a  low  level  of  pool  I. 

As  a  basis  for  the  calibration  of  the 
model,  it  was  assumed  that  in  the  data 
reported  by  Stanford  and  Smith  (1972), 
the  contribution  of  pool  I  was 
negligible,  and  the  decomposition 
constants  for  N0  were  assigned  to  pool 
II.   Values  for  the  constants  k,  h,  and  S 
(table  7.1)  were  obtained  for  pool  I  by 
fitting  the  integral  function  of  rate 
equation  1  against  the  average  carbon 
dioxide  released  from  a  variety  of 
14C-labeled  whole  cells  (Molina  et  al . 
1980),  corrected  (Paul  and  Van  Veen  1978) 
for  a  28  percent  carbon  incorporation, 
and  for  pool  II,  against  the  net 
mineralization  data  of  Stanford  (1972) 
and  Molina  et  al .  (1980).   The 
decomposition  of  sucrose  and  cellulose 
was  assumed  to  be  accounted  for  with  a 
simple  exponential  equation  (S=l) ,  and 
the  specific  rate  used  by  Beek  and  Frisel 
(1973)  was  retained. 


42 


Table  7.1. 

Kinetic  constants  of  the  substrates 


Kinetic 
constants 


Sucrose 


Substrates 


Cellulose 


Pool  I 


Pool  II 


k(day_1) 

h(day_1) 

C:N 

EFFAC 

EFFNO 

EFFHUM 


1.0 

.75 

1,000 
2.6 


1.0 

.02 

1,000 
2.6 


0.56 
.33 
.04 
^.O 
1.6 
1 .2 


0.16 
.16 
.006 
26.0 

2.6 

.0 


1  Fixed  by  calibration  of  the  model  on  the  experimental 
data  from  1st  preincubation  and  corresponding  incubation 


of  the  Chester  soil. 

2  Set  equal  to  C:N  or  EFFAC  or  pool  I. 

What  follows  is  a  brief  description  of  a 
Chichester  et  al.  (1975)  experiment.   Two 
soils,  Chester  and  Hagerstown,  were 
preincubated  with  sucrose,  cellulose 
(2,500  mg'kg"1  C  each),  and  three  levels 
of  nitrate  (50,  100,  and  200  mg'kg"1 
NO3" -1->N)  .   Five  preincubation  periods 
(approximately  1  to  5  weeks)  were 
considered.   After  each  week  of 
preincubation,  sucrose  was  added  to  the 
soil  to  balance  the  C-CO2  lost. 
Therefore,  for  each  soil,  15  different 
soil  samples  were  available,  each 
presenting  various  degrees  of  isotopic 
nitrogen  incorporation  characterized  by 
the  combination  of  5  preincubation 
periods  and  3  nitrate  levels .   Each  of 
the  30  soil  samples  was  then  incubated  to 
measure  the  potentially  mineralizable 
nitrogen  and  isotopic  nitrogen  (recently 
incorporated  nitrogen)  according  to  the 
technique  of  Stanford  and  Smith  (1972) . 


Results  from  the  first  week's 
preincubation  (carbon  dioxide  kinetics) 
and  from  the  corresponding  incubation 
(net  nitrogen  mineralization  kinetics)  of 
the  Chester  soil  for  the  three  nitrate 
levels  were  used  to  calibrate  the  C:N 
ratio  of  the  soil  active  organic  phase, 
EFFAC,  and  EFFNO   (table  7.1)  and  to 
define  the  function  /jjj  (fig.  7.4). 
Denitrif ication  data  were  used  to 
estimate  the  proportionality  constant 
(CSTDEN)  that  was  fixed  to  C  substrate 
decomposed  per  mg  of  N03~-N  that  was  not 
recovered  (25  mg  for  the  Hagerstown  soil 
and  150  mg  for  the  Chester  soil) .   The 
rate  of  nitrification  (which  did  not 
affect  the  dynamics  of  the  simulation 
since  only  denitrification  was  active 
during  the  preincubation  and  total 
inorganic  nitrogen  was  measured)  was 
arbitrarily  set  to  5  mg-kg"1  N03"-N 
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Influence  of  nitrogen  availability  on  the 
rate  of  residue  decomposition.   Function 
jijj  obtained  by  calibration  on  the 
experimental  data  from  the  first  week 
preincubation  and  corresponding 
incubation  of  the  Chester  soil. 


formed  per  day.   Nonsymbiotic  nitrogen 
fixation  was  not  used,  and  the 
corresponding  proportionality  factor  was 
not  calibrated. 

Comparison  between  the  experimental  and 
computed  data  is  shown  in  table  7 . 2  and 
figures  7.5-7.8  for  a  subjective 
validation  of  the  model  (Stanford  and 
Smith  1972) .   Computed  carbon  dioxides 
were  higher  than  the  experimental  ones 
for  the  Hagerstown  soil,  which  "owing  to 
its  fine  texture  and  poor  aggregation, 
apparently  had  a  greater  local 
accumulation  of  CO2  than  the  well -aerated 
Chester  soil"  (Stanford  and  Smith  1972) . 
In  accordance  with  the  experimental 
observations,  low  levels  (<1  mg^kg"1  N) 
of  residual  fertilizer  were  computed 
after  each  preincubation.   The  dynamics 
of  pool  I  (labile  plus  resistant)  and 
pool  II  (resistant  only)  during  the 
preincubations  are  shown  in  figures  7.9 
and  7.10. 


When  possible,  input  values  were  assigned 
as  indicated  in  Chichester's  paper.   The 

initial  concentration  of  pool  I  was 

r    1 

arbitrarily  set  at  1  mg»kg"x  C, 
equivalent  to  about  10°  bacteria  per  gram 
of  soil.   The  initial  amount  of  pool  II 


Table  7.2. 

Cumulative  amount  of  carbon  dioxide  (mg  C-CO2)  released 
per  kilogram  of  Chester  soil  during  the  initial  5  days 
of  each  preincubation 


Initial  fertilizer  nitrate  (mg'kg"1  N) 


Pre- 

50 

100 

200 

incuba- 

experi- 

simu- 

experi- 

simu- 

experi- 

simu- 

tion 

mental 

lated 

mental 

lated 

mental 

lated 

1st 

1,506 

1,428 

1,598 

1,777 

1,725 

1,700 

2d 

111 

1,212 

747 

1,075 

950 

1,125 

3d 

1,182 

900 

1,096 

843 

1,209 

967 

4  th 

1,044 

991 

953 
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Figure  7.5. 

Simulated  actual  efficiency  factor  for 
the  decomposition  of  sucrose  plus 
cellulose  and  average  production  of 
carbon  dioxide  from  200  g  of  Chester  soil 
for  the  first  and  fifth  preincubation  for 
the  200  and  50  mg^kg"1  N03"-N  initial 
additions.   Each  cross  (+)  represents  the 
experimental  average  taken  over  the 
interval  defined  by  the  cross  and  the 
preceding  one.   Simulation  averages 
(continuous  line)  are  taken  over  24-hour 
intervals .   Experimental  data  provided  by 
F.W.  Chichester. 

for  both  soils  was  estimated  to  2,650 
mg'kg--1-  C  from  an  experimental  NQ  value 
of  353,  assuming  a  C:N  ratio  of  7.5  for 
pool  II,  a  value  higher  than  the  one 
subsequently  found  by  calibration.   The 
input  tillage  flag  (TILL)  was  raised  only 
during  the  first  computational  time  step 
of  the  preincubation  and  of  the 
incubation  process ,  to  increase  the  rate 
of  mineralization  following  the  drying  up 
and  sieving  of  the  soil  sample. 
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TIME.  HOURS 


250 


Denitrif ication  was  activated  during  the 
preincubation  phases,  when  large  amounts 
of  sucrose  were  added  to  the  soil,  and 
replaced  by  nitrification  during  the 
mineralization  process.   The  amount  of 
carbon- sucrose  added  at  the  beginning  of 
each  preincubation  period  corresponded  to 
the  experimentally  determined  rather  than 
simulated  C-CO2  released  from  the 
previous  preincubation. 

Discussion 

Chichester's  et  al .  (1975)  experimental 
data  document  the  influence  of  inorganic 
nitrogen  on  the  rate  of  residue 
decomposition  and  on  the  corresponding 
efficiency  factor  for  biological 
incorporation.   Higher  initial  nitrate 
concentrations  induced  higher  levels  of 
denitrification  and  of  nitrogen 
incorporation  in  the  soil  active  organic 
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Figure   7.6. 

Computed  (continuous  line)  and 

experimental  (X)  nitrogen  lost  by 

denitrif ication  during  the  first 

preincubation. 


phase  as  later  evidenced  by  the  higher 
net  mineralization.   Yet  CO2  evolution 
during  the  same  first  preincubation 
period  was  not  correspondingly  increased. 
Sucrose  was  added  to  the  soil  at  the 
beginning  of  each  preincubation  to 
compensate  for  the  C-CO2  lost  during  the 
preceding  preincubation.   Yet  the  rate  of 
CO2  evolution  was  reduced  after  the  first 
preincubation  when  no  more  inorganic 
nitrogen  could  be  extracted  from  the 
soil.   The  structure  of  the  model  could 
not  account  for  the  CO2  and 
denitrif ication  data  until  a  variable 
efficiency  factor  and  the  reduction 
factor  /i^  were  introduced.   Attempts  to 
define  /i^j  as  a  function  of  the  soil 


5 

TIME,   (WEEK) ' 

Figure  7.7. 

Simulated  (continuous  line)  and 
experimental  (+)  cumulative  N  net 
mineralization  after  the  first,  third, 
and  fifth  preincubations  for  the  200, 
100,  and  50  mg'kg"1  15N03"-N  initial 
additions . 


inorganic  nitrogen  concentration,  as 
suggested  by  Hunt  (1977),  or  as  a 
function  of  the  rate  of  nitrogen 
mineralization  failed.   Instead,  the  use 
of  a  ratio  C  rate:N,  expressed  as  the 
potential  rate  of  carbon  residue 
decomposed  to  the  amount  of  nitrogen 
available  for  bio-incorporation,  allowed 
for  the  calibration  of  a  n^   function  with 
a  regularly  decreasing  curvature  (fig. 
7.4). 
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Figure  7.8. 

Simulated  (continuous  line)  and 

experimental  (+)  cumulative  isotopic  N 

net  mineralization  after  the  first, 

third,  and  fifth  preincubations  for  the 

200,  100,  and  50  mg.kg"1  15N03"-N  initial 

additions. 

In  addition  to  the  nitrogen  reduction 
factor,  3  out  of  the  21  internal 
parameters  (C:N  active  organic  phase, 
maximum  efficiency  factor  EFFAC ,  and 
EFFNO)  were  calibrated  on  the  first  week 
preincubation  and  corresponding 
incubation  data  from  the  Chester  soil. 
As  indicated  earlier,  pool  I  and  pool  II 
are  operational  simplifications  of  the 
morphological  and  chemical  entities  found 
in  the  active  organic  phase,  themselves 
endowed  with  a  wide  range  of  C:N  ratios 
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Figure   7.9. 

Simulated  organic   nitrogen   in  pool   I 
(total)    and  pool   II    (resistant  only) 
during  the   preincubations    for   the   200   and 
50  mg'kg"-'-    •'■5N03"-N   initial   additions. 
The   arrows    indicate  when  sucrose  was 
added  at   the  beginning  of   the   second, 
third,    fourth  and   fifth  preincubaiton 
periods.      At   time   zero,    the   soil  was 
enriched  with  sucrose   and  cellulose   at 
2,500  mg-kg"1   C  each. 

(McGill  et  al.    1981).      The   C:N  ratio  of 
the   active   organic   phase   cannot  but 
represent   an  average.       In  NCSOIL,    this 
average    is   assumed   to  be   constant   and  was 
calibrated  to   6.      Average   C:N  ratios   of 
5 . 8   and   5 . 2   corresponding   to   labeled  and 
unlabeled  hydrolyzable   C   and  N  have  been 
measured  from  a   soil   enriched  with 
acetate   and  nitrate    (McGill   et   al .    1974). 
A  variable   efficiency  factor  with  a 
maximum  of  0.6,    a  percentage    found  to  be 
satisfactory  for   the  performance   of 
NCSOIL,    is   a  well -documented  feature  of 
soil  metabolism    (Van  Veen  1977   and  McGill 
et  al.    1981).      Preliminary  runs  made  with 
a  value   of  0.1   for   the   efficiency  factor 
of  pool   II   formation  were  unsatisfactory: 
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Figure  7 . 10 . 

Simulated  organic  isotopic  nitrogen  in 
pool  I  (total)  and  pool  II  (resistant 
only)  during  the  preincubations  for  the 
200  and  50  mg'kg'l  15N03"-N  initial 
additions . 

EFFNO  had  to  be  adjusted  to  0 . 2  for  the 
model  to  give  an  optimal  account  of  the 
experiment  data.   NTRM,  in  its  present 
form,  does  not  have  a  submodel  to 
simulate  the  dynamics  of  the  soil  gas 
phase.   Therefore  variations  in  the  rate 
of  denitrification  could  not  be  expressed 
in  terms  of  changes  in  the  partial 
pressure  of  oxygen  and  carbon  dioxide. 
When  denitrification  was  activated,  the 
proportionality  constant  CSTDEN  reflected 
the  extent  to  which  the  available 
biochemical  energy  was  used  to  reduce 
nitrate .   A  much  greater  flow  of 


decomposed  carbon  was  needed  in  the 
Chester  soil  to  reduce  nitrate  than  in 
the  poorly  structured  Hagerstown  soil. 
It  is  possible  that  under  conditions  of 
water  saturation,  CSTDEN  values  lower 
than  25  may  be  reached.   As  little  as  4 
mg  C  decomposed  per  mg  NC>3"-N  denitrified 
have  been  observed  (Hilali  and  Molina 
1979). 

In  the  following,  computer  simulated 
substrate  concentrations  and  flow  rates 
were  used  to  analyze  the  soils'  behavior 
in  Chichester's  experiment.   During  the 
first  preincubation,  and  for  the  200  mg« 
kg"1  N0_3-N  treatment,  the  nitrogen 
reduction  factor  remained  high  (/i^>1.0). 
Therefore  much  energy  was  generated  from 
sucrose  catabolism  for  nitrate  reduction. 
The  combination  of  denitrification  and 
nitrate  assimilation  rapidly  reduced  the 
level  of  nitrate  in  the  soil.   After  the 
third  day  of  preincubation,  the 
concentration  was  lowered  from  200  to  10 
mg'kg"-1-  N03_-N  in  both  the  Chester  and 
the  Hagerstown  soils.   In  spite  of  this 
drastic  drop  in  nitrate  concentration, 
the  C  rate:N  ratio  was  kept  low  enough  to 
maintain  ^>1 . 0  throughout  the  first 
preincubation.   This  situation  resulted 
from  the  input  of  mineralized  nitrogen 
returned  from  a  buildup  active  organic 
phase  combined  with  a  low  rate  of  sucrose 
decomposition  along  with  rapidly 
decreasing  sucrose  concentrations.   By 
contrast,  the  50  mg»kg_1  N03~-N  treatment 
did  not  provide  enough  inorganic  nitrogen 
either  directly  as  nitrate  or  indirectly 
as  mineralized  nitrogen  to  raise  fi^   above 
its  lowest  value  0.2.   Intermediary 
values  (0.3-0.9)  were  computed  for  the 
100  mg'kg'1  N03~-N  treatment.   Therefore, 
slow  rates  of  both  immobilization  and 
denitrification  corresponded  to  low 
initial  nitrate  concentrations,  not 
solely  because  of  inorganic  nitrogen 
substrate  limitations,  but  also  because 
of  reduced  rates  of  soil  catabolic 
activity.   In  subsequent  preincubations, 
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after  all  the  initial  nitrate  had  been 
incorporated  into  organics,  the  regular 
addition  of  sucrose  was  sufficient  to 
immediately  immobilize  the  mineralized 
ammonium  nitrogen  that  was  generated  at  a 
rate  ranging  from  2.3  to  3.3  mg'kg"1  N 
per  day.   The  resulting  ammonium 
transient  state  (<5  mg'kg"1  N)  was  low 
enough  to  keep  a  high  C  rate:N  ratio, 
even  for  low  rates  of  sucrose 
decomposition  reduced  by  a  /ijj  that  did 
not  rise  above  0.2. 

Carbon  dioxide  evolution  did  not  parallel 
the  soil  catabolic  activity.   This  was 
obtained  in  the  model  by  a  variable 
efficiency  factor  that  was  calculated  to 
direct  into  pool  I  no  more  carbon  than 
needed  by  the  C:N  ratio  of  pool  I  and 
allowed  by  the  amount  of  nitrogen 
available  for  immobilization.   Thus, 
inorganic  nitrogen  that  was  released  from 
pool  I  and  pool  II  by  mineralization  and 
incorporated  in  pool  I  by  immobilization 
served  as  a  link  for  a  double  feedback 
loop  in  the  carbon  flow  between  each 
residue  component  and  pool  I. 

The  computational  time  step  was  set  at 
0.1  day  during  the  preincubation  and  1.0 
day  during  the  incubation.   However, 
subsequent  runs  showed  that  a  1 . 0  day 
time  step  could  have  been  used  as  well 
for  the  preincubation;  not  even  the 
ammonium  transient  state  was 
significantly  changed  when 
immobilization  was  solely  fed  on  a 
limited  amount  of  mineralized  nitrogen. 

Several  peculiar  features  of  the  N  net 
mineralization  experimental  data  were 
also  recorded  by  the  simulated  data  (fig. 
7.7):   No  observable  influence  of  the 
length  of  the  preincubation  on  the  net 
mineralization  kinetics,  a  longer  lag 
phase  for  lower  initial  nitrate  levels, 
more  cumulative  N  net  mineralization  from 
the  Chester  soil  than  from  the 
Hagerstown,  and  the  difference  between 
the  net  mineralization  kinetics  for  the 


50  and  100  mg.kg"1  N03"-N  initial 
addition  more  pronounced  for  the  Chester 
soil  than  for  the  Hagerstown.   At  the  end 
of  each  preincubation  period,  all  the 
sucrose  had  been  degraded;  however, 
sufficient  cellulose  was  left  to  cause 
the  lag  phase  by  immobilization  of  some 
of  the  mineralized  nitrogen.   Simulation 
runs  made  with  no  cellulose  input  at  the 
beginning  of  the  incubations  gave  a 
linear  kinetics  of  net  N  mineralization 
without  a  lag  phase  and  with  increasing 
slopes  for  increasing  initial  nitrate 
additions . 
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Chapter  8 

UNSATURATED  FLOW  SUBMODEL 

M.J.  Shaffer  and  S.C.  Gupta1 


Notation: 


6 
t 

D(0) 
K(0) 

MS 
Ksat 

^sat 
Z 

AX 

At 

i 

J 

^LB 

^UB 
FLUX 


b 

Ps 

Psl 

Pc 

*wilt 
*air 


Volumetric  water  content 

Time 

Diffusion  coefficient 

Unsaturated  hydraulic 

conductivity 

Macroscopic  sink  term  for  water 

Saturated  hydraulic 

conductivity 

Saturated  water  content 

Soil  depth 

Nodal  sparing  with  depth 

Time  step  size 

Time  indicator  (which  time  step) 

Depth  indicator  (which  node) 

Volumetric  water  content  at  lower 

boundary  (deepest  node) 

Volumetric  water  content  at  upper 

boundary  (top  node) 

Maximum  rate  of  infiltration  for 

positive  flux  or  potential 

evaporation  rate  for  negative 

flux 

Soil  matric  potential  (suction) 

Air  entry  potential 

Campbell  "b"  coefficient;  slope 

of  a  log-log  plot  of  i>(B)   vs.  6 

Percent  sand 

Percent  silt 

Percent  clay 

Volumetric  water  content  at 

permanent  wilt 

Volumetric  water  content  at  ipr 


Introduction 

The  subject  of  unsaturated  flow  of  water 
in  soils  has  received  considerable  and 
detailed  study  over  the  past  four 

1  Shaffer  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108,  and 
Gupta  is  an  associate 
professor,  Soil  Science 
Department,  University  of 
Minnesota,  St.  Paul,  MN  55108. 


decades.   Various  models  have  been 
proposed  to  simulate  the  infiltration  and 
redistribution  processes.   These  models 
range  from  simple  piston  displacement 
models  (Dutt  and  Tanj i  1962),  through 
semiempirical  approaches  such  as  the 
Green  and  Ampt  equation  (Green  and  Ampt 
1911,  Li  et  al.  1976),  to  analytical  and 
numerical  solutions  of  the  Richard's 
equation  (Warrick  et  al .  1971,  Hanks  and 
Bower  1962) . 

The  unsaturated  flow  submodel  reported 
here  is  a  numerical  solution  to  the 
Richard's  equation  (Dutt  et  al .  1972). 
Simultaneous  infiltration,  evaporation, 
redistribution,  deep  percolation,  and 
plant  uptake  of  water  are  simulated  as  a 
function  of  time.   Upward  and  downward 
movement  of  water  are  allowed. 

Submodel  inputs  include  precipitation 
and/or  irrigation  amounts  and  dates , 
potential  or  maximum  rates  of 
infiltration  and  evaporation,  potential 
transpiration,  and  current  values  for 
crop  root  distribution,  soil  water 
contents,  and  soil  hydraulic  properties. 
Model  outputs  include  soil  water  contents 
and  fluxes  with  time  and  depth,  soil 
evaporation  and  infiltration  rates ,  the 
transpiration  rate,  and  the  rate  of  deep 
percolation  (may  be  positive  or 
negative) .   The  submodel  also  outputs 
values  for  the  RATI01  coefficient,  which 
is  the  ratio  of  actual  to  potential 
transpiration.   This  is  used  to  regulate 
crop  growth  in  water  stress  situations. 


Model  Description 

This  submodel  uses  a  finite  difference 
solution  to  the  moisture  flow  equation  in 
one  dimension  to  simulate  the  waterflow 
processes  discussed  previously.   The 
basic  flow  equation,  which  includes  a 
macroscopic  sink  term  for  plant  water 
uptake  (MS) ,  is  written, 

30/at=a|D(0)30/3z-K(0)]/3z-MS.       [1] 
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The  finite  difference  solution  to 
equation  1  can  be  written  as 

3      3  L  j+l/2v  j+1  j+1  j   j 


--^U-^i^K1 


,i-l 


+2AXK 


i-1 
J -1/2 


]/2AX 


2  MS1. 
J 


J-l   J'* 

[2] 


The  water  flux  between  adjacent  nodes  is 
calculated  for  each  time  step,  using  the 
relationship 


FLUX= 


Tk1" 


1/2     j+l/2v   j+1     j+1 


.1-1, 


/2AX] 


At. 


Values  for  D(0)  are  calculated  using  an 
estimate  of  the  water  content  (6)    at  the 
end  of  the  next  time  step  and  averaged 
over  two  depth  nodes .   This  approximation 
for  anticipated  water  content  is  written 


j+l/2    j+1  j+1  j   j   " 


+0' 


[4: 


The  upper  boundary  (soil  surface)  is 
defined  in  terms  of  a  water  content  (#ub) 
that  varies  over  time,  depending  on  the 
surface  conditions.   The  values  used  for 
#UB  in  the  finite  difference  equations 
are  determined  using 

0UB=4=(2.OAX/D^/2).(Hm-Kji;1/2) 


+1.8^"1+^"1: 


[6; 


where 

FLUX  =  the  maximum  rate  of  infiltration 
for  positive  water  flux  (water  on 
the  surface)  or 


FLUX  =  the  potential  evaporation 

rate  for  the  case  where  water  is 
evaporating  from  the  soil  surface. 

(Note:   the  assumption  is  made  that 
evaporation  is  zero  when  infiltration  is 
taking  place) . 


The  actual  values  computed  for  FLUX  (+  or 
- )  depend  on  the  soil  water  contents  and 
hydraulic  properties.   This  allows  the 
soil  infiltration  rate  to  vary  with  time 
during  an  infiltration  event  and  the  soil 
evaporation  rate  to  vary  as  the  soil 
dries  out. 


where  Y  is  a  weighting  factor  defined  as 
Y=0.7(Ati)/Ati'1.  [5] 

Boundary  Conditions 

The  lower  boundary  is  defined  in  terms  of 
a  constant  water  content  (#lb) .   The 
value  used  for  fl^B  can  either  be  the 
saturated  theta  value  (0 s)  or  some  other 
value  less  than  saturation.   The  value 
used  must  be  determined  based  on  the  best 
approximation  of  conditions  in  a 
particular  field  situation. 


Soil  Water  Characteristic  Curves 

The  values  for  soil  hydraulic 
conductivity  (K)  and  diffusivity  (D)  as  a 
function  of  soil  water  content  (8)   must 
be  provided  either  with  a  set  of  user 
supplied  equations  or  by  using  the 
Campbell  method  for  the  soil  water 
characteristic  curves  (Campbell  1974) . 
This  method  assumes  the  moisture  release 
curve  (suction  vs.  theta)  can  be 
represented  by  the  equation 


iK0)«iM0/0sat) 


7] 
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and  the  hydraulic  conductivity  vs.  theta 
relationship  by 


K(0)=KsatWsat) 


2b+3 


[s; 


Diffusivity  (D(0))  is  defined  by  the 
relationship 

D(0)=K(0)/O0/dV), 


which  gives 

D(0)«(btf  -K    •9h+2)/dh+l 
r   sat        sat 


[9] 


10 


when  equations  7,  8,  and  9  are  combined. 

The  relationship  given  by  equation  7 
plots  as  a  straight  line  on  log-log 
paper.   The  slope  is  equal  to  minus  the 
empirical  constant  b, 


b=-ln(V>2/V>i)/ln(02/0l) 


11] 


If  the  user  does  not  have  data  for  rj>   vs . 
0,    the  model  will  calculate  a  value  for  b 
based  on  the  soil  particle  size  analysis, 
bulk  density,  and  percent  organic  matter 
using  a  set  of  regression  equations 
developed  by  Rawls  et  al .  (1981).   The 
same  equations  are  used  to  estimate  new 
soil  water  characteristics  curves 
whenever  the  above  properties  change. 
This  provides  one  of  the  links  between 
tillage  and  basic  soil  properties. 

In  addition  to  a  value (s)  for  the 
empirical  constant  b,  the  user  must 
provide  values  for  air  entry  potential 
(V»r)  and  saturated  hydraulic  conductivity 
(Ksat) .  r/>v   can  be  derived  from 
experiment  or  computed  from  an  empirical 
equation  based  on  the  particle  size 
analysis.   In  the  latter  case,  the  model 
will  use  the  equation 


A  value  for  Ksaj-  must  be  provided  by  the 
user.   This  should  be  based  on 
information  obtained  under  field 
conditions  so  that  at  least  one  point  on 
the  hydraulic  conductivity  vs.  water 
content  curve  is  field-verified. 

The  model  also  requires  values  for  6 sat 
and  the  water  contents  at  permanent  wilt 

(^wilt)  anc*  alr  ^rv  (^air)  •   Users  can 
supply  these  values  from  their  own  data 
or  allow  the  program  to  estimate  them 
based  on  the  regression  equations. 


Model  Sensitivity  Analyses 

Figures  8.1-8.3  contain  plots  of  soil 
water  content  (0)  vs.  depth  and  time. 
These  graphs  present  typical  results 
obtained  from  the  unsaturated  flow 
submodel.   They  illustrate  the  types  of 
responses  the  submodel  makes  to  water 
applications  and  evapotranspiration. 

In  figure  8.1,  a  fixed  water  table  was 
maintained  at  the  lower  boundary. 
Several  deep  percolation  events  are 
visible  that  contributed  water  to  the 
aquifer.   The  high  variability  of  0   near 
the  soil  surface  is  characteristic  of 
many  soils.   The  infiltration  and 
subsequent  evapotranspiration  and 
redistribution  of  precipitation  events 
causes  the  patterns  shown. 

Figures  8.2  and  8.3  are  depth  plots  of 
soil  water  content  at  different  times 
after  an  infiltration  event.   Note  the 
rapid  downward  movement  of  the  soil- 
wetting  front  in  figure  8.2.   The  effects 
of  unsaturated  flow  after  the  front  has 
passed  are  also  visible  at  166.4  days. 


V>r=0.147/(10PS+3.3PSL+Pc).34,000    [12]    Model  Verification 


to  estimate  V*r  if  a  value  is  not 
provided. 


The  finite  difference  solution  to  the 
unsaturated  flow  equation  has  seen 
extensive  use  and  verification  over  the 
years  by  a  number  of  researchers  and 


53 


i&'°    i\f€  ( 


Figure  8.1. 

Theta  vs.  depth  over  time. 
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Figure  8.2. 

Moisture  profile  changes  as  12.8  cm  of 

water  is  applied  at  day  166.0. 
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Figure  8.3. 

Soil  moisture  on  Julian  dates  143  and 

174. 


others  (Dutt  et  al.  1972,  Nimah  and  Hanks 
1973,  Wierenga  et  al .  1974).   The 
particular  solution  technique  used  here 
has  been  tested  by  the  U.S.  Bureau  of 
Reclamation,  the  U.S.  Environmental 
Protection  Agency,  the  U.S.  Department  of 
Agriculture,  and  several  universities. 
Examples  of  the  type  of  agreement 
obtained  between  predicted  output  and 
observed  data  are  given  in  figures  8.4- 
8.6.   The  model  gives  good  agreement  with 
experiments  for  predicted  water  contents 
with  depth  (and  time) ,  infiltration  into 
the  soil,  and  deep  percolation  leaving 
the  bottom  boundary.   In  addition,  use  of 
model  predictions  for  water  fluxes  and 
water  contents  over  depth  and  time  have 
been  shown  to  adequately  predict  the 
movement  and  distribution  of  solutes 
(Wierenga  et  al.  1974). 
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Figure  8.4. 

Observed  and  predicted  soil  moisture 

tension  for  lysimeter  2  at  37-  and  73-cm 

depth. 
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Figure    8.5. 

Predicted  and  observed  cumulative 
infiltration  with  time. 
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Figure  8.6. 

Observed  and  calculated  cumulative 

drainage  effluent  curves. 
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Chapter  9 

CROP  GROWTH  SUBMODEL  (CORN) 

M.J.  Shaffer1 


Introduction 

Simulation  models  for  the  growth  of 
agricultural  crops  are  being  rapidly 
developed  by  several  researchers.   Most 
of  these  models  contain  computer  code 
that  is  specific  for  each  individual 
crop,  although  some  contain  general 
relationships  common  to  many  crops .   In 
general,  software  availability  dictates 
what  crop -specific  top  growth  submodels 
can  be  used  in  simulation  models  of  the 
crop -soil -water  continuum.   A  more 
general  crop  model  would  be  highly 
desirable  for  use  in  large  simulators. 
We  are  currently  developing  a  multicrop 
submodel  for  use  in  NTRM.   This  software 
will  allow  the  simulation  of  about  20 
agricultural  crops. 

For  initial  NTRM  modeling  purposes,  we 
selected  the  Nebraska  corn  growth  model, 
CORNGRO  (Childs  et  al .  1977,  Tscheschke 
and  Gilley  1980) .   This  submodel  is 
somewhat  specific  for  corn  but  does 
contain  some  relationships  that  are 
adaptable  to  other  crops.   The  corn  crop 
was  selected  because  of  its  agronomic 
importance  in  the  Midwest.   The  CORNGRO 
model  was  kept  intact  with  respect  to 
ideal  growing  conditions  as  a  function  of 
solar  radiation.   Interactions  with  the 
environment  were  modified  or  added  to 
simulate  the  effects  of  air  and  soil 
temperature,  soil  water,  soil  aeration, 
soil  strength,  soil  nitrogen,  and  soil 
salinity.   The  resulting  submodel  for 
corn  growth  was  then  interfaced  with  the 
overall  NTRM  model.   Key  linkages  include 
the  uptake  of  soil  NH4+-N  and  N03~-N,  the 
DEFNIT,  RATIO,  and  RATI01  parameters,  and 
other  input  parameters  such  as  air 
temperature  and  solar  radiation. 

1  Research  soil  scientist, 
U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 


Model  Description 

The  corn  growth  submodel  operates  as  a 
producer  and  consumer  of  carbohydrates. 
Carbohydrates  are  produced  by 
photosynthesis,  which  is  simulated  in  the 
model  by  the  relationship, 


P=PSNK» FLEAF* FRAD. FCHO'RATIOl 


[I! 


where 
P 
PSNK 


net  photosynthesis  (mg 
CH20.dm"2»h"1), 


gross  photosynthesis  at 
1.4  ly/min 
=  65.0  mg  CH20'dm-2.h-1, 
FLEAF  =  fraction  of  total  leaf  area 
that  is  sunlit  (0-1) , 

FRAD  =  2.36R/(1+1.65R) ;  R  is  solar 

radiation  in  ly/min  (this  is  an 
adjustment  function  for  gross 
photosynthesis  based  on 
incoming  solar  radiation)  (0-1), 

FCHO  =  (84.5-A)/84.5;  A  is  available 
carbohydrate  (this  factor 
sets  maximum  available 
carbohydrate  at  84.5  mg  CH2O 
•cm"2)  (0-1),  and 
RATI01  =  actual  transpiration/ 

potential  transpiration 
(0-1). 

Carbohydrates  are  consumed  in  the  model 
by  the  respiration  processes  for  growth 
and  maintenance.   Maintenance  respiration 
is  calculated  using  the  expression 


RM=C(0.044+1.9.10-3-T+1.10-3.T2)TDMP 


[2] 


where 

RM  = 

C  = 

T  = 
TDMP  ■ 


maintenance  respiration 
(mg  CH20'dm"2'h"1) , 
maintenance  coefficient 
(0.74  mg  CH20«dm"2«h"1«g"1) , 
air  temperature  (°C),  and 
total  dry  matter  produced  to 
date  (g  dry  matter/plant) . 
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Growth  respiration  is  computed  by 
calculating  total  respiration  and  then 
subtracting  maintenance  respiration.   For 
total  respiration, 

RTOT=(3.857»10-3+3.829'10"3.T)A.DEFNIT 
•RATIO,  [3] 

where 

RTOT  =  total  respiration  (mg  C^Odm"2 
•h"1), 
DEFNIT  =  actual  N  uptake/N  uptake  demand, 
and 
RATIO  =  root  mass  produced/potential 
root  mass  production; 


for  growth  respiration, 

RG=RT0T-RM, 
where 


:*] 


RG  =  growth  respiration  (mg  C^O'dm"2 


•h"1) 


Crop  growth  is  computed  from  growth 
respiration  using 

DMP=3.76RG, 

where 


:s] 


DMP  =  growth  produced  (mg  C^O'dm 
•h_1). 


-2 


Available  carbohydrate  (A)  is  updated  on 
an  hourly  (1/20  day  in  NTRM)  basis  using 


A=At.1+P- RTOT -DMP, 

where 

At.]_  =  available  carbohydrate  from  the 
previous  time  step. 


6] 


Total  dry  matter  production  in  the  tops 
is  computed  using  the  expression 

tasseling 
DRYMAT=X3.58RG.LA»FRACT/1,000  (g  dry 
emergence 

matter/plant),  [7] 

where 

o 

LA  =  leaf  area  (dm~^) ,  and 

FRACT  =  fraction  of  dry  matter  in  the 
tops . 

Total  dry  matter  production  in  the  ears 
is  given  by 

maturity 
DRYEAR=£3.58RG»LA'FRACT/1,000  (g  dry 
tasseling 

matter/plant).  [8] 

Leaf  area  from  emergence  to  tasseling  is 
computed  daily  using 

LA=LA0CD ( 1+AK1 • RGDA»  FRACT • FLFDRY) ,   [ 9 ] 

where 

o 
LAOCD  =  previous  day's  leaf  area  (dm^) , 

AK1  =  AK2-ATEMP(35-T)/25.AK2, 

RGDA  =  daily  sum  of  growth  respiration, 

FLFDRY  =  ratio  of  leaf  growth  to  total 

dry  matter  production  (varies 

from  0.9  at  emergence  to  0.0  at 

tasseling) , 

AK2  =  crop  coefficient  (determined 

from  calibration) ,  and 

ATEMP  =  temperature  sensitivity 

coefficient  (0.72). 
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After  tasseling,  leaf  area  is  computed 
using  the  equation 

LA=AALA(1- (DDAY- (DDT+100) )/(DDM- (DDT 
+100)),  [10] 

where 

AALA  =  leaf  area  at  tasseling, 
DDAY  =  the  number  of  accumulated  degree 
days , 
DDT  =  number  of  degree  days  to 

tasseling,  and 
DDM  =  number  of  degree  days  to  maturity, 

Note  that  in  equation  10,  the  value  for 
AALA  is  constant  and  100  growing  degree 
days  are  assumed  to  occur  before  a 
significant  reduction  in  LA  begins. 

The  conversion  from  dry  matter  produced 
(DRYMAT  or  DRYEAR)  to  uptake  of  soil 
nitrogen  (NH4+-N  +  N03~-N)  is  done 
according  to  the  following  scheme: 


A  general  flowchart  for  the  crop  growth 
submodel  is  presented  in  figure  9.1. 
Note  that  each  of  the  major  growth 
processes  has  its  own  subroutines.   These 
are  called  from  the  DRYMAT  subroutine, 
which  in  turn  is  called  from  the 
administrative  routine  CROPGRO.   Input 
data  for  the  submodel  are  assembled  by 
the  CROPG  routine. 

Examples  of  typical  submodel  output  are 
shown  in  figures  9.2-9.4.   Note  that  the 
submodel  begins  the  corn  growth  simula- 
tion when  the  crop  has  emerged  from  the 
soil  and  the  leaf  area  is  1  dm  .   Leaf 
area  increase  continues  until  tasseling, 
where  leaf  senescense  and  ear  development 
begins.   These  two  processes  continue 
until  crop  maturity.   At  this  point,  the 
crop  growth  submodel  is  no  longer  called 
and  any  crop  materials  remaining  on  the 
surface  are  treated  as  surface  residues. 


TNAPTK=NFRACT • DRYM 


Davs   after 

N   fraction 

emergence 

drv  matter 

0-15 

0.0221 

16-30 

.0173 

31-45 

.0150 

46 -maturity 

.0120 

Accumulated  degree  days  are  used  to 
estimate  the  dates  of  tasseling  and 
maturity.   The  following  equation  is  used 
to  make  the  degree  day  calculations  based 
on  hourly  (or  l/20th  day)  air 
temperatures : 


DDAY=£(T-10)/20. 

The  summation  is  limited  to  air 
temperatures  in  the  range 

10°C<T<35°C. 


11] 
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Figure  9.1. 

Crop  growth  submodel  (corn) . 
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Typical  Output  and  Model  Sensitivity 

Figures  9.2,  9.3,  and  9.4  contain  plots 
of  typical  dry  matter,  grain,  and  leaf 
area  output  data,  respectively,  from  the 
crop  growth  submodel.   Time  zero  is  just 
after  emergence  with  a  leaf  area  of  1 
dm  .   Note  that  grain  development  (fig. 
9.3)  begins  at  tasseling  and  continues 
until  maturity.   Grain  production  before 
maturity  represents  the  cumulative  dry 
matter  to  date  incorporated  into  the 
developing  ear.   Leaf  area  development, 
(fig.  9.4)  progresses  until  tasseling, 
where  leaf  senescense  begins . 


Model  Verification 
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Verification  of  the  corn  growth  submodel 
has  been  done  by  Childs  et  al .  (1977), 
Tscheschke  and  Gilley  (1980),  and  by  me 
using  the  CORNGRO  model  as  modified  for 
use  in  the  NTRM  model.   Figures  9.5  and 
9.6  contain  typical  results  obtained  by 
Childs  and  Tscheschke  using  their 
respective  versions  of  the  CORNGRO  model 
The  observed  data  are  from  experimental 
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Figure   9.3. 

Grain  production   (corn) 
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Leaf  area  vs.    time    (corn). 
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Figure   9.2. 
Dry  matter  vs.    time    (corn). 
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Figure   9.5. 

Actual  and  predicted  leaf  area  and  dry 
matter  production  for  Sandhills 
Agricultural  Laboratory,  Nebraska,  1975. 


plots  in  Nebraska.   Note  that  the  overall 
approach  generated  reasonably  good 
agreement  with  experiment  for  Nebraska. 

The  CORNGRO  model  version  used  in  the 
NTRM  model  was  tested  using  data  from 
four  sites  in  the  Midwestern  United 
States,  including  two  sites  in  Minnesota 
(Apple  Valley  and  Lamberton) .   A 
comparison  of  predicted  and  observed  corn 
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Figure   9.6. 

Yield  comparisons   on  CORNGRO  model    for 
various    irrigation   treatments   and 
varieties   at  Mead  and  Sandhills 
Agricultural   Laboratory,    Nebraska. 


yields   obtained   from   this   analysis    is 
presented   in   figure    9.7.      The   predicted 
results   were   obtained   from  runs   using   the 
overall   NTRM  model  with   submodel 
interactions   considered.      The   scatter   in 
the   data  points    is   probably  typical   of 
results    that   can  be   expected   from   the 
NTRM  and   similar  models. 
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NTRM  model- -corn  yields. 
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Chapter  10 

ROOT  GROWTH  SUBMODEL 

M.J.  Shaffer  and  C.E.  Clapp1 


Introduction 

A  one -dimensional  dynamic  submodel  for 
root  growth  was  developed  to  simulate 
root  penetration,  root  branching,  and 
root  death.   The  submodel  was  designed  to 
serve  as  an  interactive  component  of  the 
NTRM  model.   Initially,  model  coeffi- 
cients were  developed  for  corn,  but 
enough  flexibility  was  retained  in  the 
model  to  allow  the  introduction  of 
coefficients  for  different  crops. 
Modeled  environmental  effects  on  root 
growth  include  soil  temperature,  soil 
bulk  density,  soil  aeration,  soil 
strength,  soil  water  potential,  soil 
salinity,  and  soil  nitrogen.   The  model 
simulates  root  branching,  death,  and 
penetration  rates  as  a  function  of  these 
variables  and  if  necessary,  redistributes 
the  root  growth  within  the  rooting  zone 
to  exploit  the  most  favorable  soil 
conditions . 

Baseline  or  reference  rates  of  root 
growth  (branching  and  penetration)  and 
root  death  are  taken  as  functions  of  crop 
growth  stage,  and  dry  matter  production 
by  the  plant.   Because  plant  dry  matter 
production  is  in  part  regulated  by 
environmental  conditions  encountered  by 
the  roots,  the  tops  and  roots  are 
interconnected  by  a  continuous  feedback 
loop. 

Root  density  changes  are  the  result  of 
root  growth  and  root  death  (slough) . 
Root  slough  becomes  input  into  the  NTRM 
nitrogen  and  carbon  transformation 
submodel.   Root  death  rates  are  a 
function  of  growth  stage,  soil 
environmental  conditions,  and  the 
availability  of  root  mass. 


1  Research  soil  scientists, 
U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 


Physical  and  chemical  barriers  (or 
partial  barriers)  may  or  may  not  alter 
the  total  length  (mass)  of  roots 
produced.   In  these  situations,  the  model 
attempts  to  alter  the  root  growth 
pattern,  within  limits,  and  grow  the 
roots  in  a  more  favorable  part(s)  of  the 
root  zone  (Taylor  1974) .   The  model  is 
constrained  by  maximum  root  densities, 
maximum  upper  limits  on  growth  rates,  and 
adverse  environmental  conditions  else- 
where in  the  root  zone.   When  the  model 
cannot  redistribute  the  root  mass,  a  loss 
of  roots  occurs  and  a  signal  is  sent  to 
the  top  growth  submodel  via  the  RATIO 
coefficient  to  reduce  dry  matter 
production.   This  is  analogous  to  a 
hormonal  feedback  mechanism  in  real 
plants- -when  cytokinin  production  by  the 
roots  is  inhibited  by  stress  (Skene 
1975). 

Indirect  impacts  on  crop  growth  are  also 
possible.   Any  change  from  the  ideal  case 
will  alter  the  root  distribution  with 
depth  and  time.   This,  in  turn,  may  have 
indirect  impacts  on  root  and  top  growth 
due  to  changes  in  water  and  nutrient 
uptake  patterns . 

A  general  flowchart  for  the  root  growth 
submodel  is  shown  in  figure  10.1.   The 
ROOTS  subroutine  calls  the  computational 
subroutines  and  calculates  the  root 
slough  and  root  distribution  for  each 
soil  segment.   BDTEMH  sets  up  the 
environmental  variables  for  each  1-cm 
root  horizon  (fig.  10.2),  while  the  EXTN 
and  PSEUD  routines  are  involved  with  root 
redistribution  and  extension  in  stress 
cases . 

Root  Branching  Submodel  (RDEMOD) 

This  subroutine  computes  the  rate  of  root 
branching  as  a  function  of  the 
environmental  variables.   It  calls 
function  routines  for  soil  temperature 
(T) ,  soil  strength  (S) ,  soil  aeration 
(A),  soil  potential  (includes  salinity) 
(SS) ,  and  soil  nitrogen  (N) . 
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Figure  10.1. 

Root  growth  submodel . 

The  rate  of  root  branching  as  a  function 
of  these  variables  can  be  written 

t-br  sbr 

Rbr=Rbr0-J(aRbr/3T)dT-J(3Rbr/5s)ds 


tbj 


>br, 


SS 


br 


Abr 
-J(3Rbr/3A)dA-f(3Rbr/3SS)dSS 

Abr0         ssbr0 

Nbr 
-J(3Rbr/3N)dN-IA  terms. 

Nbr, 


1] 


Each  function  subroutine  evaluates  its 
respective,  integral  term  in  equation  1. 
The  branching  rate  given  by  Rbr  is  the 
rate  under  reference  conditions  that  are 
taken  as  being  ideal  with  respect  to  root 
growth. 


"I -DIMENSIONAL 
ROOT  GROWTH 


1  cm2 


SOIL  SURFACE 


II         |l 


RATE  OF 
EXTENSION 


1  cm 


RATE  OF  DENSITY 
INCREASE  OR  DECREASE 


150  cm  MAXIMUM 
DEPTH  OF  ROOTING 

Figure   10.2. 

Root  growth  model  configuration. 

Equation  1  separates  the  environmental 
effects  on  root  growth  into  a  series  of 
independent  and  cross-product  IA 
(interaction)  terms.   This  approach 
differs  from  the  methods  used  by  some 
modelers  of  selecting  the  process  with 
the  greatest  negative  impact  on  rate  as 
being  rate  determining,  or  representing 
the  rate  as  a  product  of  a  reference  rate 
and  a  series  of  multiplicative 
environmental  variables  normalized 
between  0  and  1.   From  mathematical  and 
theoretical  standpoints,  equation  1  is 
the  most  comprehensive  and  rigorous  of 
the  three  methods.   Evaluation  of  the 
interactive  terms  presents  some 
difficulties  that  can  be  solved  using  a 
combination  of  field  research  and 
modeling. 
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The  root  growth  submodel  allows  the  use 
of  any  one  of  the  above  methods  by 
setting  an  appropriate  control  switch. 

To  simplify  the  calculations  and  relate 
the  three  modeling  methods,  Rbr  can  be 
written  as  a  function  of  the  reference 
branching  rate  (Rbr  )  anc*  a  factor  (f) 
that  varies  between  zero  and  one. 


Rbr_fmin* Rbr0  - 


[6] 


Rbr=f,Rbr0- 

By  rewriting  equation  1,  we  have 

fcbr 
Rbr=Rbr0+Rbr0'J(3f/3T)dT 

Cbr0 
sbr  Abr 

+Rbr  -fof/as)ds+Rbr  .f(af/aA)dA 


2] 


Sbi 


Ab] 


SSbr  Nbr 

+Rbro»J(af/ass)dss+Rbro.J(af/aN)dN, 


ss 


br. 


Nbr, 


3] 


where  the  integral  J*(af/3Vari)dVari  for 
each  variable  (Vari)  is  equal  to 

d-fVari) 
Simplifying,  we  have, 

N 
Rbr=Rbr0-Rbr0'X(l-fvari).  and      W 

N 
f=l-X(l-fVari)  [5] 

1 

for  the  general  case  where  a  series  of 
independent  and  cross-product  terms  is 
being  used. 

For  the  modeling  method  that  sets  f  equal 
to  the  minimum  fyari  >    ecluation  2  becomes 


and  for  the  multiplicative  case  we  have 


Rbr=n.fVari.Rbr 


Root  Extension  Submodel  (REXMOD) 


[7] 


Root  extension  (penetration)  is  simulated 
using  a  modified  form  of  equation  3.   The 
terms  for  soil  matric  potential  (includes 
salinity)  and  soil  nitrogen  were  omitted, 
and  different  values  were  used  for  the 
reference  rates  (see  table  10.1).   The 
assumptions  were  made  that  root 
penetration  rates  are  not  a  function  of 
soil  nitrogen  levels  or  soil  matric  plus 
osmotic  potentials.   Experimental 
evidence  suggests  that  roots  will  readily 
penetrate  into  dry  (or  saline)  soils  and 
soils  devoid  of  nitrogen,  provided  water 
and  N  are  available  elsewhere  in  the  root 
zone.   Note,  however,  that  the  soil 
strength,  soil  temperature,  and  soil 
aeration  terms  were  retained. 

Root  Growth  Factor 

The  factor  for  soil  temperature  is 
calculated  using  the  relationship  for 
corn  shown  in  figure  10.3. 

The  response  to  changes  in  temperature  is 
assumed  to  be  linear  from  15 °C  up  to 
25°C.   A  plateau  is  then  taken  from  25°C 
to  30 °C  where  a  linear  reduction  is  made 
to  0  at  40°C  (Allmaras  et  al .  1964, 
Taylor  et  al.  1972,  Nielson  1974). 

The  submodel  is  designed  to  generate 
compensatory  root  growth  in  response  to 
temperatures  that  are  not  in  the  optimal 
range  (Russell  1977)  .   This  growth  is 
restricted  only  by  a  maximum  rate  of 
growth  and  an  upper  limit  on  root  density 
both  of  which  can  be  set  by  the  user. 

The  function  used  for  soil  strength, 

f (strength) ,  contains  the  variables  bulk 
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Table   10.1. 

Root  coefficients   for  corn,    Rosemount,    MN    (1981) 


Days  from 

Zero 

order 

1st  order 

emergence 

(cm 

/day) 

(day" 

-1) 

K 

K1 

K2 

K3 

(extension) 

(branching) 

(b: 

ranching) 

(death) 

0-35 

1.54 

0.01 

0.05 

0.02 

35-50 

8.00 

.01 

.05 

.02 

50-80 

4.00 

.01 

.05 

.02 

80-90 

.00 

.00 

.01 

.02 

90+ 

.00 

.00 

.00 

.02 

1  Valid  for  root  densities  <0.1  cm/cnr 

2  Valid  for  root  densities  >0.1  cm/cnr 

3  Literature   estimates. 


IDEAL  TEMP 


(25°-30°C) 


Figure   10.3. 

Soil   temperature   function. 

density    (BD)    and   soil  matric   potential 
(V>)    and   is   written 

f  (strength)  — 1.095+0.  9735 •  10" 7V>2 

-  0.1084*  10  "2V>+4.636BD 

-2.044BD2  [8] 

(BD   range=l.  1-2.0,    tJ>  range=0-700   cm). 


A   graphical   display   of   this    function   is 
shown   in   figure   10.4.      The  basic   data 
were   adapted   from  Taylor   and  Gardner 
(1963) ,    and  the  working  functions  were 
derived  using  a  second  order  regression 
approach. 

The  value   for  BD  used  in  equation  8    is 
obtained   from 


BDL=(SS-90.0)'KP+BDLo,    and 


BD=BD0+BDL0-BDL,  or 


BD=BDo-(SS-90.0).KP, 
where 


[9] 
10] 

11] 


BDQ  =  the  soil  bulk  density, 
KP  =  the  penetration  coefficient  (range 

0.004-0.020), 
SS  =  the  sum  of  percent  sand  and 
percent  silt,  and 
BDL0  =  the  limiting  bulk  density  for  the 
reference  soil. 

These  equations  incorporate  the  work  of 
O'Connell  (1975)  to  extend  Taylor  and 
Gardner's  data  for  a  fine  sandy  loam  soil 
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Figure   10.4. 

Soil  strength  function. 

to  other  soil  textures.   In  equation  9, 
the  KP  value  is  the  slope  of  the  equation 
relating  bulk  density  to  SS  (percent  of 
particles  >0.05^).   The  90.0  and  BDLQ 
refer  to  the  SS  and  limiting  bulk 
densities,  respectively,  for  Gardner's 
soil.   The  BDLQ  value  can  be  computed  as 
a  function  of  suction  by  using  equation  8 
and  setting  f( strength)  equal  to  0. 
However,  it  cancels  out  of  the  final 
expression  for  BDQ  equation  11.    In 
effect,  equations  9,  10,  and  11  shift  the 
x-axis  location  in  figure  10.4  to  account 
for  changes  in  soil  texture. 

The  soil  aeration  function,  f (aeration) , 
is  shown  graphically  in  figure  10.5. 
This  function  is  set  equal  to  unity  (1.0) 
in  the  soil  water  content  (range  0  to  85 
percent)  of  saturation.   A  linear 
reduction  is  then  taken  between  this 
point  and  saturation  where  the  function 
becomes  equal  to  0. 


The  function  for  soil  salinity,  f(salin- 
ity) ,  is  shown  in  figure  10.6.   This 
relationship  looks  at  the  sum  of  the 
matric  and  osmotic  potentials.   Soil 
water  potential  values  for  10,  25,  and  50 
percent  growth  reductions  were  adapted 


from  Bernstein  (1964)  and  U.S.  Department 
of  Agriculture's  Salinity  Laboratory 
Staff  (1954).   The  assumption  was  made 
that  the  function  equals  zero  at  50  atm. 
The  salinity  function  alters  the  rates  of 
root  branching,  and  the  submodel  attempts 
to  increase  root  density  in  more 
favorable  portions  of  the  root  zone.   The 
root  extension  rates  are  not  affected  by 
this  function.   The  function  for  soil 
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Figure  10.5. 

Soil  aeration  function. 
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nitrogen,  f (nitrogen) ,  is  shown  in  figure 
10.7.   This  equation  generates  an 
increased  rate  of  root  branching  in  those 
portions  of  the  root  zone  containing  the 
most  ideal  concentrations  of  NH^-N+NC^" 
-N  (Maizlish  et  al .  1980).   The  function 
also  effectively  shuts  off  the  rates  of 
root  branching  and  extension  where  toxic 
levels  of  NH4  -N  are  encountered  (Pesek 
et  al.  1971). 

In  all  cases,  the  root  growth  will  be 
redistributed  to  the  most  desirable 
portions  of  the  soil  root  zone.   If  this 
cannot  be  done ,  the  hormonal  feedback 
mechanism  is  activated  to  restrict 
overall  dry  matter  production. 

Compensatory  Root  Growth 


POTENTIAL 

COMPENSATORY  GROWTH 
(PCG) 


DEPTH 
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Figure  10.8. 

Compensatory  root  distribution. 


Root  growth  is  redistributed  by  the  model 
in  response  to  adverse  conditions  with 
respect  to  soil  temperature,  soil 
strength,  soil  water  potential,  soil 
aeration,  and  soil  nitrogen  levels.   The 
algorithm  used  to  accomplish  this  is 
illustrated  in  figure  10.8.   Note  that 
the  total  potential  loss  in  root  growth 
for  each  1- cm- depth  element  is  -  computed 
for  each  daily  time  step.   The  ability  of 
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the  active  root  zone  to  compensate  for 
this  loss  is  calculated  for  each  element 
and  summed  over  the  root  zone.   This  sum 
is  compared  with  the  total  potential  loss 
to  determine  the  total  compensatory  root 
growth  and  the  total  actual  loss.   The 
actual  compensatory  growth  is  then 
proportioned  over  the  depth  elements 
according  to  their  ability  to  grow  the 
additional  roots.   Any  remaining  growth 
loss  is  used  in  the  computation  of  the 
RATIO  coefficient. 

The  compensatory  growth  potential  is 
calculated  as  a  user- supplied  fraction  of 
the  maximum  growth  rate . 


Root  Distribution  With  Depth 

The  fraction  distribution  (Fr)  of  the 
roots  with  time  and  depth  is  computed 
using  the  following  relationship  for  each 
daily  time  step: 


N 
Fr^=Length^/£  Lengthy. 
i=l 


12 


Figure  10.7. 

Soil  nitrogen  function. 
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The  root  length  summation  is  done  from 
the  soil  surface  down  to  the  current 
rooting  depth.   This  gives  Fr^  the 
definition  of  current  fraction  of  total 
root  length  contained  in  depth  element  i . 

Relationship  of  Root  Weight  and  Root 
Length 

To  compute  the  root  weights  needed  for 
comparison  with  top  weights  and  total  dry 
matter  production  and  to  use  in  computing 
the  root  slough  weights  needed  by  the 
nitrogen  and  carbon  transformation 
submodel,  an  expression  was  derived  based 
on  corn  root  data  collected  at  Rosemount, 
MN.   All  tillage,  residue,  and  nitrogen 
treatments  from  the  1981  growing  season 
were  combined  to  produce  a  single  data 
set  for  root  length  versus  weight.   The 
assumption  was  made  that  the  subsurface 
portion  of  the  plant  stem  was  not  part  of 
the  root  system.   The  root  core  samples 
were  weighted  according  to  the  volume  of 
soil  they  represented.   A  best-fit 
multiple  linear  regression  equation  was 
developed  from  these  data  as  follows : 


Wgt=2 . 076+8 . 925Lth-8 . 47 8 J  Lth 


13] 


where 

Wgt  =  the  root  weight  (mg/cnr)  and 
Lth  =  the  root  length  (cm/cm-*)  . 

The  R^  value  for  the  regression  was  0.826 
with  153  cases  analyzed. 

Equation  13  was  developed  using  root 
lengths  in  the  range  0.175  to  1.71 
cm/cm-* .   Root  length  values  less  than 
0.175  should  be  used  in  the  following 
equation: 


Wgt=0.520Lth 


14 


This  represents  a  straight  line  from  0  up 
to  the  0.175  length.   Equation  14  was 
developed  to  avoid  computing  negative 


numbers  for  root  weight  in  this  region, 
given  the  form  of  the  best- fit  regression 
(equation  13) . 

A  combined  plot  of  predicted  versus 
observed  root  weights  obtained  from  the 
regression  equation  and  equation  14  is 
shown  in  figure  10 . 9 .   The  R^  value  for  a 
linear  1:1  correlation  is  0.847. 

Model  Calibration  and  Sensitivity 

The  root  growth  submodel  was  calibrated 
using  1980  and  1981  corn  root  data  from 
the  USDA-ARS  research  plots  at  Rosemount, 
MN.   Based  on  these  data  and  data  from 
the  literature,  rate  coefficients  were 
developed  for  root  extension,  root 
branching,  and  root  death  for  various 
stages  of  corn  development.   The  values 
for  these  coefficients  are  summarized  in 
table  10.1. 

Observed  and  calculated  values  for  root 
density  on  several  dates  during  the  1980 
growing  season  are  shown  in  figure  10.10. 
These  data  represent  the  moldboard  plow 
tillage  treatment.   Note  that  the 
observed  root  densities  are  reproduced 
with  reasonably  good  accuracy. 

Figures  10.11  and  10.12  illustrate  a 
model  sensitivity  comparison  in  three 
dimensions  for  somewhat  typical  root 
growth  pattern  (figure  10.11)  and  a 
situation  where  a  partially  restrictive 
pan  is  present  (figure  10.12).   Note  the 
striking  differences  in  root  growth  over 
time  for  these  two  cases. 

Changes  in  root  distribution  over  a 
growing  season  are  illustrated  in  figure 
10.13.   Each  line  in  the  figure  starting 
from  top  to  bottom  represents  the 
fraction  of  total  roots  in  each 
succeeding  15 -cm- depth  increment  down  to 
75  cm. 
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Figure  10.9. 

Predicted  versus  observed  root  weights. 
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Observed  and  calculated  root  densities, 

Rosemount,    MN. 
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Figure   10.11. 

Typical  root  length  vs.  depth  over  time, 
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Typical  root  length  vs.  depth  over  time 
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Root  growth  as  a  function  of  depth  and 

time  (corn). 
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Chapter  11 

TILLAGE  AND  SURFACE-RESIDUE-SENSITIVE 

POTENTIAL  EVAPORATION  SUBMODEL 


D.R.  Linden,  D.M.  Van  Doren,  Jr. , 
J.K.  Radke,  and  D.C.  Reicosky1 

Introduction 

Evaporation  is  the  process  by  which  water 
is  returned  to  the  atmosphere.   Evapora- 
tion from  the  soil  is  of  particular 
interest  to  this  model,  although  other 
aspects  of  evaporation  will  also  be 
considered.   It  is  affected  by  the  energy 
available  to  heat  and  vaporize  water,  the 
ease  with  which  the  vapor  can  move  away 
from  the  soil,  and  the  ease  with  which 
water  will  move  to  the  evaporating 
surface  from  within  the  soil.   Evapora- 
tion from  a  freshly  wetted  soil  is 
generally  limited  by  the  available  energy 
and  aerodynamic  transfer  conditions 
(energy  limiting  or  constant -rate  phase) 
until  the  surface  soil  becomes  dry  enough 
that  waterflow  within  the  soil  will  not 
meet  the  atmospheric  demands  (soil 
limiting  or  falling-rate  phase). 
Describing  or  predicting  evaporation  thus 
involves  describing  waterflow  within  the 
soil,  with  special  care  being  given  to 
the  boundary  condition  at  the  soil 
surface.   Tillage  affects  both  soil 
waterflow  properties  and  the  surface 
boundary  conditions  and  thus  evaporation. 


1  Linden  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108;  Van 
Doren  is  a  professor  emeritis, 
Ohio  Agricultural  Research  and 
Development  Center,  Wooster, 
OH  44691;  Radke  is  a  research 
soil  scientist,  U.S. 
Department  of  Agriculture, 
Agricultural  Research  Service, 
Rodale  Research  Center, 
Kutztown,  PA  19530;  and 
Reicosky  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Morris,  MN 
56267. 


Tillage  is  a  disturbance  of  soil  to  some 
depth.   The  depth  of  tillage  depends  on 
the  objective  of  the  disturbance. 
Shallow  tillage  is  used  to  control  weeds, 
to  prepare  seedbeds,  to  place  seeds 
within  the  soil,  to  mix  an  additive  (for 
example,  herbicides  or  fertilizers)  with 
the  soil,  to  change  the  shape  of  the 
soil -atmosphere  boundary,  to  incorporate 
residues,  or  to  change  the  heat  and  water 
transport  properties.   Deep  tillage  may 
be  used  for  any  of  the  above  reasons  and 
to  break  up  a  root-restricting  or  water- 
impeding  layer  within  the  soil.   This 
paper  describes  how  these  tillage 
disturbances  affect  potential  evaporation 
rate.   The  effects  of  closely  associated 
phenomena,  such  as  crop-residue 
condition,  are  also  considered. 

Potential  evaporation  has  been  separated 
and  isolated  from  evapotranspiration 
because  tillage  affects  the  surface 
boundary  and  the  soil  properties  and  thus 
directly  affects  evaporation,  while  in 
many  ways  tillage  affects  transpiration 
more  indirectly,  for  example  by  changing 
the  depth  distribution  of  plant  roots, 
crop  stand,  or  weed  populations.   The 
effects  of  tillage  on  the  soil  waterflow 
process  have  been  presented  elsewhere 
(Linden  1982)  and  will  not  be  discussed 
here . 

Tillage  has  been  shown  to  increase  short- 
term  evaporation  rates  (Holmes  et  al. 
1960,  Allmaras  et  al .  1977)  and  to 
decrease  longer  term  evaporation  (Willis 
and  Bond  1971,  Gill  et  al.  1977). 
Surface  residues  are  known  to  decrease 
short-term  evaporation  (Van  Doren  and 
Allmaras  1978) .   Although  the  observed 
effects  of  tillage  are  somewhat 
contradictory,  they  are  consistent  with 
some  popular  concepts  about  the  use  of 
tillage  to  control  evaporation.   In  humid 
regions,  tillage  is  considered  a 
mechanism  for  enhancing  evaporation  and 
thus  speeding  soil  drying.   In  more  arid 
regions,  shallow  surface  tillage  is 
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Aerodynamic  transfer  coefficient 
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considered  a  mechanism  for  water 
conservation.   The  contradictory  nature 
of  these  concepts  clearly  establishes  the 
need  for  predictive  capability.   Soil, 
climate,  the  timing  of  tillage 
operations,  the  type  of  tillage  tool,  and 
the  depth  of  tillage  are  interactive 
factors  that  affect  the  soil  water 
system.   Surface  residues  are  generally 
considered  as  a  mechanism  for  conserving 
water.   A  model  of  evaporation 
considering  short-  and  long-term  effects 
can  help  in  making  tillage  decisions  on  a 
universal  basis,  which  would  be 
especially  useful  where  the  balance 
between  short-  and  long-term  effects 
varies  with  the  season. 

Methods  of  predicting  potential 
evaporation  range  from  largely 
statistical  approaches  to  theoretical 
descriptions  of  soil  water  and 
atmospheric  flow  systems  (van  Bavel  1966, 
Gardner  1973,  van  Bavel  and  Hillel  1976). 
The  statistical  approaches  cannot 
generally  be  adapted  for  variations  in 
soil  condition  and  are  thus  not  useful 
for  the  purposes  of  this  paper.   It  is 
not  our  intent  in  this  paper  to  review, 
discuss,  and  judge  the  merits  of  various 
existing  approaches  but  rather  to  present 
some  new  concepts  on  how  any  model  might 
be  adapted  to  tillage  and  residue 
conditions.   The  model  includes  these 
effects  of  tillage:   (1)  soil-atmosphere 
boundary  effects  on  evaporation  energy 
resulting  from  radiation  trapping  and  (2) 
residue  effects  on  energy  and  water 
transfer  in  the  residue  mulch  layer. 

The  model  is  intended  to  be  used  as  an 
evaporation  model  throughout  the 
nonfrozen  soil  year  and  thus  includes 
growing  crop  canopy  effects  on:   (1) 
light  transmission  to  the  soil  surface 
and  (2)  the  aerodynamic  transfer  of  water 
vapor  through  the  canopy.   A  conceptual 
diagram  of  the  processes  and  effects 


considered  in  the  evaporation  submodel  is 
shown  in  figure  11.1.   The  concepts 
represented  in  the  diagram  will  be 
discussed  in  the  following  pages. 
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Figure  11.1. 

Conceptual  diagram  of  the  processes  and 
functions  in  the  tillage- res idue- 
sensitive  evaporation  submodel. 
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Potential  or  Energy-Limiting  Evaporation 
Flux 

Energy- limiting  evaporation  rate  is 
described  by  a  combination  equation  as 
presented  by  Penman  et  al.  (1967) 
including  the  Businger  (1956)  transfer 
function  expressed  as 


pe=[a/7(rn-h-g+rNl)+v.dt-c] 

/[A/7+l], 


[i: 


where 
A  = 

7  = 

RN  = 

H  = 

G  = 

RNL  = 

V  = 
DT  = 

C  = 

PE  = 


the  slope  of  the  vapor  pressure  vs 

temperature  curve , 

the  psychrometric  constant, 

the  net  shortwave  radiation, 

the  energy  used  to  heat  air, 

the  energy  used  to  heat  soil, 

the  net  longwave  radiation, 

the  vapor  deficit, 

the  L  aerodynamic  transfer 

coefficient  for  water  vapor, 

a  combination  of  constants  (see 

below) ,  and 

the  potential  evaporation  rate . 


Equation  1  is  used  to  calculate  the 
potential  evaporation  rate  from  the 
measured  input  data  of  total  solar 
radiation,  maximum  and  minimum  air 
temperature,  and  total  daily  wind  travel 
with  intermediate  equations  relating  the 
measured  input  data  to  the  terms  of 
equation  1.   The  intermediate 
relationships  will  be  discussed  below. 
Equation  1  is  used  to  calculate  the 
potential  evaporation  flux  across  three 
planes  of  reference:   (1)  the  top  of  the 
crop  canopy,  (2)  the  top  of  the  residue 
mulch,  and  (3)  the  soil  surface.   The  net 
radiation,  Rjj,  and  aerodynamic  transfer 
coefficient,  D-j,  terms  of  equation  1  are 
modified  to  include  the  effects  of 
tillage,  residues,  and  canopy  at  each  of 
these  planes. 


Net  Radiation 

Net  radiation,  Rjj,  may  be  measured 
directly  or  calculated  from  measured 
total  radiation  and  the  albedo  of  the 
evaporating  surface.   An  alternative 
input  could  be  pan  evaporation  corrected 
for  wind  effects  and  converted  to  an 
"effective  radiation  term." 

Albedo  of  the  surface  includes  the 
effects  of  soil  roughness  and  water 
content  and  the  crop -residue  mulch 
reflectivity  and  water  content. 

The  potential  evaporation  rate  is  very 
sensitive  to  the  net  radiation  available 
for  evaporation  and  thus  to  the  albedo  of 
the  surface  (Idso  et  al.  1975,  van  Bavel 
and  Hillel  1976) .   Soil  roughness  reduces 
the  reflectivity  of  a  soil  (Bowers  and 
Hanks  1965,  Allmaras  et  al.  1972,  Cary 
and  Evans  1975,  Gausman  et  al .  1977).   In 
1978,  Linden  presented  the  concept  and  a 
model  to  account  for  the  reduced 
reflectivity  by  the  mechanism  of 
radiation  trapping  on  roughened  soil 
surfaces  (Linden  1978,  Linden  1979,  Cruse 
et  al.  1980). 

Briefly,  the  model  assumes  that  reflected 
sunlight  or  emitted  terrestrial  radiation 
is  perfectly  scattered  and  that 
variations  in  topography  (roughness)  can, 
therefore,  intercept  this  radiation  so 
that  it  does  not  escape  from  the  soil 
system  into  the  atmosphere.   Microtopo- 
graphic  data  and  reflectivity  for  a 
smooth  condition  can  be  used  to  estimate 
the  magnitude  of  radiation  trapping  and 
reflectivity  on  the  rough  soil.   The 
model  was  tested  by  comparing  measured 
and  predicted  reflectivity  for  a  4-hour 
midday  period  on  a  rough  soil  when  the 
soil  surface  was  dry.   Reflectivity  of  a 
smooth  tilled  soil  was  measured  to  be 
0.179  for  use  in  predicting  the 
reflectivity  of  the  rough  surface . 
Reflectivity  for  the  rough  surface  was 
obtained  from  an  upright  and  an  inverted 
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Epply  pyranometer.   The  predicted 
reflectivity  was  0.131,  and  the  measured 
was  0.130. 

Since  the  model  appeared  to  give 
reasonable  predictions  of  reflectivity,  a 
general  relationship  between  roughness 
and  daily  average  reflectivity  was 
developed. 

Several  sets  of  microtopography  data  with 
various  roughness  conditions  were 
subjected  to  model  calculations,  and  the 
predicted  reflectivities  were  correlated 
with  the  standard  deviation  of  the 
microtopographic  data.   The  standard 
deviation  of  elevations,  s,  is  used  as  an 
index  of  the  roughness  with  higher  s 
occurring  on  rougher  surfaces.   Predicted 
daily  average  reflectivities  relative  to 
a  smooth  condition  at  a  similar  water 
content  plotted  as  a  function  of 
roughness  index  are  shown  in  figure  11.2. 
As  roughness  increases,  reflectivity 
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Figure   11.2. 

Predicted  relative  albedo  and  potential 
evaporation  rate  (PE)  as  a  function  of 
roughness  index  (standard  deviation  for 
elevation  measurements) .   Data  shown  is 
relative  to  values  for  smooth  surface 
(roughness  index  =  0). 


decreases  (fig.  11.2),  which  will  tend  to 
increase  the  potential  evaporation  rate, 
PE.   For  purposes  of  illustrating  this 
effect  of  tillage  on  evaporation, 
predicted  PE's  relative  to  a  smooth 
condition  with  roughness -affected  albedo 
as  the  only  variable  input  are  also  shown 
in  figure  11.2.   The  PE's  were  estimated 
with  equation  1  using  reasonable  and 
nonvariant  assumptions  for  the  remaining 
equation  terms.   Roughness -related  albedo 
did  not  have  a  large  impact  on  the 
potential  evaporation  rate  (fig.  11.2). 
However,  small  differences  in  this  rate 
can  have  a  large  impact  on  cumulative 
evaporation  over  long  periods. 

Reflectivity  of  the  soil  is  also  affected 
by  the  water  content  of  the  surface . 
Drier  soils  will  reflect  more  radiation 
than  moist  soils.   The  effect  of  changing 
volume  and  therefore  area  (along  the 
surface)  fractions  of  water  and  air  on 
the  reflectivity  of  the  soil  surface  can 
be  expressed  as  a  simple  weighted  mean 
between  water  and  soil  components  and 
written  as 


aw=(0.06r+(0.5-r)ao)/0.5, 
where 


[2: 


a^j  =   the  reflectivity  of  the  soil, 
r  =  the  volumetric  water  content  of  the 

surface ,  and 
aQ  =  the  reflection  coefficient  of  dry 

soil. 

As  r  increases,  a  higher  proportion  of 
the  surface  is  occupied  by  water,  which 
has  a  reflection  coefficient  of  about 
0.06.   The  reflectivity  of  soil,  aw,  has 
limits  of 


0.55ao<aw<ao 


[3: 


The  reflectivity  of  the  soil,  a^,  and  the 
roughness  index,  s,  are  used  as  inputs  to 
determine  the  albedo  of  the  soil  surface 
as  affected  by  radiation  trapping 
discussed  previously  but  that  cannot  be 
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reduced  to  a  simple  functional  equation 
because  of  the  effect  of  the  altitude 
angle  of  the  sun  on  radiation  trapping 
(Linden  1979) .   An  effective  daily  albedo 
of  the  soil  is  thus  a  function  of  time  of 
year  and  is  recomputed  daily  by 
considering  the  altitude  angle  of  the  sun 
in  the  evaporation  submodel.   The 
reflectivity  of  a  dry  smooth  soil,  aQ,  is 
a  necessary  input  parameter  to  the  model. 

The  reflectivity  of  residues  is  affected 
by  age  and  water  content.   The  effect  of 
age  of  the  residues  was  calculated  with 
an  equation  similar  to  one  given  by  von 
Hoyningen-Huene  (1971) .   The  reflecti- 
vity of  the  residue  mulch,  aR,  is  given 
by  " 


aR=ao(1.06+aRl/ao-1.06)e-°-0255t), 


[4; 


where 

ao  = 

t  = 

aRi  = 


the  reflectivity  of  dry  smooth 
soil , 

the  residue  age  in  days,  and 
the  reflectivity  of  the  fresh 
residue  mulch. 


As  time  (age  of  residues)  increases,  the 
residues  approach  the  reflectivity  of  the 
soil  material,  whereas  fresh  residues 
often  have  a  reflectivity  twice  as  high 
as  soil  material.   This  reduction  in 
reflectivity  is  due  to  both  decomposition 
of  the  residue  material  and  the  coating 
of  residue  pieces  by  splashed  soil 
material.   The  effect  of  residue  water 
content  on  the  reflectivity  of  residues, 
aR,  is  modeled  as  reducing  to  75  percent 
of  the  aR  whenever  the  residues  contain 
evaporable  water  (see  later  section  on 
water  retention  and  evaporation  from  crop 
residues) .   A  combined  surface  albedo  was 
computed  by  considering  the  fraction  of 
the  soil  surface  covered  and  not  covered 
by  residue  material.   The  fraction  of 
exposed  soil  surface,  B,  was  calculated 
from  a  refitted  exponential  equation  of 


data  presented  by  Sloneker  and 
Moldenhauer  (1977)  and  given  as 

B=e-^, 
where 


[5] 


M  =  the  dry  weight  of  residues  per  unit 

of  area  and 
f  =  a  term  dependent  on  type  of 

residues  present. 

An  overall  reflection  coefficient,  as , 
was  calculated  from  the  relative 
proportion  of  the  surfaces  covered  with 
residues,  1-B,  and  uncovered,  B,  and  the 
respective  reflectivities.   The  equation 
is  given  as 


as=aR(l-B)+awB. 


6] 


On  a  daily  basis,  the  residue  coverage  is 
recomputed  from  updated  surface  residue 
amounts  supplied  to  this  model  from  the 
nitrogen  transformation  submodel  (Chapter 
7,  "Carbon  and  Nitrogen  Transformations 
in  Soil,  Submodel  NCSOIL").   The  residue 
reflection  characteristics  and  the  soil 
reflectivity  are  recomputed.   Potential 
evaporation  thus  has  dynamic  reflection 
characteristics  that  depend  on  water 
content  of  the  soil,  roughness,  residue 
age,  and  residue  water  content.   The 
reflection  coefficients  of  fresh  dry 
residues  and  smooth  dry  soil  surfaces  are 
the  input  reflection  characteristics  from 
which  the  dynamic  reflection  properties 
are  determined. 

The  total  radiant  energy  available  for 
evaporating  water  from  the  residue  or 
soil  surface  is  computed  by  considering 
the  light  transmitted  through  a  crop 
canopy  to  the  soil  surface.   A  light 
transmission  fraction,  F^ ,  is  computed 
from  Ritchie  (1972) : 


FT=e 


-0.6L 


for  L<Ll,  or 


[7; 


F    -0.384LL.e-0.216L  for  ^       [g] 


where 

LL=0.404.1nS+1.49,  [9] 

where 

S  =  the  ground  area  occupied  per  plant 

and 
L  =  the  leaf  area  index. 

Equations  7 ,  8 ,  and  9  are  used  to 
calculate  F-p  from  leaf  area  index  and 
plant -spacing  data.   The  light 
transmission  fraction,  F^,  times  the 
total  incoming  solar  radiation  gives  the 
radiation  reaching  the  residue  and  soil 
surface . 


A  potential  evaporation  equation  (eq.  1) 
has  been  made  partially  sensitive  to 
tillage  and  residues  by  considering  the 
effects  on  albedos  and  thus  net  radiation 
available  for  evaporation.   Total 
incoming  radiation  is  the  primary 
variable  but  is  not  always  available; 
therefore  an  attempt  was  made  to  compute 
effective  radiation  from  pan  evaporation 
data,  which  is  more  readily  available. 
The  radiation-evaporation  equation  (eq. 
1)  is  rearranged  to  give  effective  radi- 
ation, R(Rjg-H-G+Rjg  )  ,  and  is  written  as 


R=[pEp(l+W/q)-V.DTc]/(W/q)0.93(l-a), 


13] 


where 


An  overall  canopy  albedo,  a,  is  computed 
from  the  weighted  radiation  fraction  as: 


a=asFT+ac(l-FT), 

where 

as  -  the  surface  albedo , 
ac  =  the  crop  albedo ,  and 
a  =  the  overall  field  albedo. 

The  net  radiation,  RN  ,  available  for 
evaporation  at  the  top  of  the  canopy 
would  be  given  as 


10 


RNc=(l-a)R, 
and  the  net  radiation,  R 


N< 


[11] 
available  for 


evaporation  at  the  soil  (or  residue) 
surface  would  be  given  as 


PEp  =  pan  evaporation, 

a  =  the  combined  albedo  of  the  field, 
0.93  =  a  crop  coefficient  factor, 

and  the  remaining  variables  are  as 
previously  defined  (eq.  1). 

The  effect  -of  canopy  interception  of 
radiation  in  reducing  the  net  radiation 
at  the  soil  surface  and  the  potential 
evaporation  rate  is  shown  in  figure  11.3. 
Net  radiation  was  calculated  for  a  range 
of  L  values  from  equation  12  and  used  in 
equation  1  as  the  only  variate  to 
develop  the  curve  in  figure  11.3.   As 
leaf  area  index  increases ,  PE  decreases 
at  a  declining  rate.   Modeled  potential 
evaporation  is  very  sensitive  to  the 
shading  caused  by  the  canopy,  as  should 
be  expected. 


RNs=d-as)FTR 


[12] 


where 


R  = 


the  total  incoming  solar 
radiation, 
the  net  radiation, 
a  and  as  =  the  weighted  mean  albedos,  and 
F<p  =  the  light  transmission 
fraction. 


RN  = 


Soil  and  Air  Heating 

The  energy  consumed  in  heating  soil,  H, 
and  air,  G,  is  assumed  to  be  0  since  no 
data  were  found  to  quantify  these 
relationships.   The  assumption  of  0 
sensible  heatflow  to  air  and  soil  is  an 
approximation  that  is  used  extensively  in 
estimating  potential  evaporation  rate. 
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Figure  11.3. 

Predicted  relative  potential  evaporation 
rate  PE  and  the  fraction  radiation 
transmitted  to  the  surface  FT  as  a 
function  of  leaf  area  index. 


Net  Longwave  Radiation 

The  net  longwave  radiation,  Rjj  ,  is  given 
by  an  empirical  calibration  equation  to 
daily  maximum  and  minimum  air 
temperatures  by  Ritchie  (1972)  and  is 
written  as 


RNl=( -2. 61+0. 00414. TL1-8)R/RC 
where 

TL=((Tmax+Tmin)/2)-7.75, 


14 


15 


Tmax  =  the  maximum  daily  air  temperature, 
Tmin  =  t^e  minimum  daily  air  temperature, 
R  =  the  measured  incoming  solar 

radiation  and  Rc  is  the  maximum 
possible  incoming  solar  radiation. 

Net  longwave  radiation  is  generally  a 
small  term  on  a  daily  average  basis. 


Delta 

The  slope  of  the  vapor  pressure  vs. 
temperature  curve,  A,  is  determined  from 
the  average  of  the  maximum  and  minimum 
temperature  for  the  day  by  an  equation 
presented  by  Ritchie  (1972): 

A_^(21. 255-5, 30VTk)].5i304/Tk2i       [16] 


where 


Tk=(Tmax+Tmin/2)+273. 


[17 


Gamma 


The  psychrometric  constant,  7,  is 
determined  from  the  temperature  of 
condensation,  which  is  assumed  to  be  at 
the  minimum  air  temperature  for  the  day 
by  an  empirical  equation  taken  from 
Ritchie  (1972).   The  ambient  air  pressure 
is  corrected  for  elevation  before 
determining  the  psychrometric  constant 


7=6.6.10-4.P(l+.00115.Tnn-T1) 


mm> 


18 


where 


Lmm 


the  elevation- corrected  average 

atmospheric  pressure  and 

the  minimum  daily  air  temperature 


Vapor  Deficit 

The  vapor  deficit,  V,  was  calculated  from 
the  algorithm  given  by  Richards  (1971)  by 
assuming  that  the  soil  surface  was  at  the 
assumed  temperature  of  condensation 
(minimum  air  temperature)  and  the  ambient 
air  vapor  pressure  could  be  calculated 
from  the  daily  mean  air  temperature  (TMAX 
+TMIN)/2.   Vapor  pressures  were  obtained 
from  the  Richards  (1971)  algorithm.   The 
vapor  pressure  deficit  was  calculated  as 
the  difference  in  the  vapor  pressure  at 
these  two  temperatures. 
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Aerodynamic  Transfer  Coefficient 

The  aerodynamic  transfer  coefficient,  D-p, 
was  calculated  from  the  daily  wind  travel 
by  an  equation  relating  the  wind  travel; 
the  logarithm  ratio  of  the  height 
measurement,  k^',    an  equivalent  roughness 
height,  Z0;  and  a  length  of  path  for 
diffusion  transfer,  Zp,  to  aerodynamic 
transfer.   The  aerodynamic  transfer 
coefficient  is  thus  calculated  from  daily 
wind  travel  using 

DT=(U.C)/[(ln(ZA+Zc+Z0/Z0+ZD))2 

+(0.1681.U.ZD/DA)]  [19] 

where 

Opek2/P=1.168.1(r3.0.622.0.412/l,013 


-1.2212-10"7, 


[20] 


Zo  = 

U  = 


the  height  of  wind  measurement, 
the  diffusion  coefficient  for  water 
vapor  in  still  air, 
a  roughness  parameter, 
the  daily  wind  run, 


and  Zc ,  ZQ ,  and  Zp  are  discussed  more 
fully  below. 

The  aerodynamic  transfer  of  water  vapor 
is  made  sensitive  to  tillage,  residue, 
and  crop  canopy  conditions  by  adjustments 
to  the  Zc  and  Zp  terms  of  equation  19  and 
not  the  ZQ  term  as  discussed  by  Lettou 
(1969) .   The  height  adjustment  factor  for 
the  wind  profile  in  the  canopy,  Zc,  was 
calculated  as 


ZC=0.6HC 


21] 


where  Hc  is  the  height  of  the  canopy. 
Equation  21  assumes  that  the  wind  as 
measured  at  ZA  height  could  have  been 
measured  at  0.6«HC+ZA  height  over  the 
canopy.   Furthermore,  the  wind  profile 


assumed  to  be  identical  so  that  the  wind 
speed  is  0  at  0.4»HC  height  and  that 
vapor  transfer  below  this  height  is  by 
diffusion  only  so  that 


ZD=0.4HC. 


[22 


The  effective  diffusion  path  length,  Zp 
is  increased  to  include  the  effects  of 
the  residue  mulch  when  considering 
potential  evaporation  from  the  soil 
surface . 


ZD=0.4Hc+(d.A(B+A)/l+(A/4)) 
where 

A=0.0127732M/dpR  and 
B=e'A 


23 


[24; 


25 


where 

d  = 

A  = 

B  = 


PR 


the  diameter  of  residue  pieces, 

the  total  vertical  cross  section 

of  residues, 

the  fraction  of  the  surface  not 

covered  by  residues, 

the  residue  dry  weight,  and 

the  density  of  residues. 


is 


Equations  21  through  25  are  used  to  make 
adjustments  in  the  aerodynamic  transfer 
coefficient,  D-p  (eq.  18),  to  include  the 
effects  of  the  canopy  and  the  residue 
mulch  on  the  evaporation  of  water.   A 
summary  of  the  radiation,  R,  albedo,  a, 
and  height  adjustment  terms  Zc  and  Zp 
used  for  calculating  potential 
evaporation  at  three  planes  is  given  in 
table  11.1. 

The  modeled  effect  of  crop  height  as  the 
only  variate  on  potential  evaporation  is 
shown  in  figure  11.4.   The  effect  of  crop 
height  was  modeled  with  equations  21  and 
22  for  Zc  and  Zp  values  used  in  equation 
19  to  calculate  D-p,  which  in  turn  is  used 
in  equation  1  to  calculate  PE.   Crop 
height  has  most  of  its  effect  on  PE  in 
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Table  11.1. 

Summary  of  parameters1  used  in  equations  11,  12,  and  19 
to  calculate  the  net  radiation  available  for  evaporation 
and  the  aerodynamic  transfer  coefficient 


Potential 

evaporation 

plane 


Radiation   Albedo     Effective  Diffusion 

reaching  height  path 

plane  of  wind  length 

measurement 


Top  of 
canopy 


R 


Top  of 
mulch 


FTR 


ZA+0 . 6HC     0 . 4HC 


Soil 
surface 


FTR 


ZA+0.6HC     0.4Hc+d.A(B+A)/ 
l+(A/4) 


1  See  text  for  an  explanation  of  variables 

the  first  meter  of  growth  (fig.  11.4)  and 
thereafter  has  little  effect.   The  effect 
of  residue  mulches  on  potential 
evaporation  for  various  environmental 
conditions  is  shown  in  figure  11.5.   The 
effects  of  residue  reflectivity  and 
diffusive  flow  resistance  are  included  in 
the  modeled  PE  rates  of  figure  11.5. 
Residues  decrease  potential  evaporation 
at  declining  rates  as  the  amount  of 
residues  increase.   The  effects  of 
residues  are  most  pronounced  at  high 
windspeeds . 


Interface  of  Potential  Evaporation  Rate 
to  Soil  Waterflow  Model 

The  potential  evaporation  algorithm  (eq, 
1)  is  operated  on  a  daily  basis  with 
modified  albedos  and  aerodynamic  height 
terms  as  the  management -affected 
variables.   The  model  is  driven  with 
radiation,  maximum  and  minimum  air 
temperatures ,  and  wind  run  as  the 
climatic  inputs.   Pan  evaporation  data 
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Figure  11.4. 

Predicted  relative  potential  evaporation 

rate  as  a  function  of  crop  height. 
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Figure   11.5. 

Predicted  relative  potential  evaporation 
rates  as  a  function  of  residue  amounts 
and  surface  residue  coverage. 

can  be  used  as  substitute  input  because 
they  are  convertible  to  an  effective 
radiation  amount  for  use  in  equation  1. 
Effective  sunshine  percentage  may  also  be 
used  as  alternative  input  since  the  model 
includes  an  algorithm  for  calculating  the 
radiation  flux  density,  which  can  be  used 
as  a  substitute  for  measured  radiation. 
The  daily  output  PE  of  this  model  is  the 
control  for  surface  water  flux  until  the 
soil  surface  can  no  longer  meet  this 
demand.   After  this,  the  actual  soil 
evaporation  rate  is  a  function  of  the 
soil  flow  properties  and  is  calculated  in 
the  waterflow  model. 

Evaporation  From  the  Crop  Canopy  and 
Residue  Mulch 

Both  living  and  dead  plant  material  above 
the  soil  surface  will  intercept  rainfall 
that  is  not  available  for  infiltration 
into  the  soil.   This  retained  water  is 
also  subject  to  evaporation.   Evaporation 


from  the  canopy  is  considered  to  occur  at 
a  rate  equal  to  the  potential  flux  rate 
at  the  top  of  the  canopy  minus  the  flux 
rate  at  the  surface ,  which  is  at  the 
bottom  of  the  canopy.   Likewise,  the 
evaporation  rate  from  the  residue  mulch 
is  considered  to  occur  at  a  rate  equal  to 
the  difference  in  potential  rate  at  the 
top  and  bottom  (soil)  of  the  mulch  layer. 
Water  is  added  to  the  canopy  and  mulch 
during  rainfall  and  is  discussed  in 
Chapter  14,  "Simulation  of  Interception, 
Surface  Roughness,  Depression  Storage, 
and  Soil  Settling." 


Summary  and  Conclusions 

In  summary,  the  model  presented  here 
would  require  the  following  as  the 
minimum  input  information,  in  addition  to 
some  climatic  data:   roughness  index; 
reflectivity  of  the  soil  in  a  smooth,  dry 
condition;  and  reflectivity  of  the 
residues.   The  output  of  the  model 
includes  the  evaporation  rate, 
evaporation  amount,  soil  water  storage, 
and  water  contents  at  specific  times. 

The  major  benefit  from  modeling 
evaporation  has  been  the  isolation  and 
study  of  single  variables  affected  by 
tillage  and  residues.   These  isolation 
studies  have  shown  that  roughness  is  of 
some  significance,  but  that  the  hydraulic 
properties  of  the  soil  are  of  major 
concern  (Linden  1982) .   Other  effects 
still  need  to  be  considered  for  a  model 
to  be  completely  sensitive  to  tillage  and 
residues.   The  effect  of  tillage -induced 
roughness  on  the  aerodynamic  transfer  of 
water  vapor  away  from  the  soil  surface; 
the  effect  of  roughness  on  nonuniform 
drying  (van  Bavel  and  Hillel  1976)  caused 
by  nonuniform  radiation  inputs;  the 
effect  of  roughness  and  porosity  on 
thermal  gradients  that  induce  waterflow; 
and  the  effects  of  porosity  on  water 
vapor  movement  through  diffusion,  mass 
air  movement,  and  soil  water  hysteresis 
should  be  included. 
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Chapter  12 
TRANSPIRATION  SUBMODEL 

M.J.  Shaffer1 


Introduction 

The  TRANSP  submodel  computes  potential 
transpiration  based  on  pan  evaporation, 
leaf  area,  and  root  distribution.   The 
potential  values  at  each  node  are  used  by 
the  unsaturated  flow  model  in  its 
macroscopic  sink  term.   The  actual 
predicted  values  for  transpiration  at 
each  node  depend  on  the  soil  hydraulic 
properties,  the  soil  water  contents,  and 
the  potential  transpiration. 

The  expression  used  to  calculate 
potential  transpiration  is  written 


PT^RTK  •  PAN .  •  ALEAF .  - 
111  1 


COEF 


[i: 


where 


PT  =  potential  transpiration  for  the 

Jth  layer, 
RD  =  the  fraction  of  roots  in  the  Jth 
layer, 
PAN  =  the  pan  evaporation  for  the  ith 
time  step, 
ALEAF  =  the  leaf  area  for  the  ith  time 
step,  and 
COEF  =  a  crop  coefficient,  which  depends 
in  part  on  the  stand  density. 

Equation  1  takes  into  account  the  total 
energy  available  for  evaporation, 
transpiration,  and  the  size  of  the  plant. 
Both  the  root  distribution  and  the  leaf 
area  are  dynamic  properties  that  vary 
with  time  according  to  what  is  going  on 
in  the  overall  system. 

The  model  extracts  water  down  to  the 
permanent  wilting  point  (usually  15  bars 
suction)  before  water  stress  begins  and 
actual  predicted  transpiration  becomes 
less  than  the  potential  predicted  value. 


1  Research  soil  scientist, 
U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 


The  model  currently  attempts  to 
redistribute  the  water  uptake  pattern 
when  stress  occurs  in  only  part  of  the 
root  zone.   Unsaturated  flow  of  water  can 
also  occur,  from  wet  to  dry  regions. 

When  water  stress  is  occurring,  a 
coefficient  called  RATI01  becomes  equal 
to  a  number  in  the  range  1>RATI01>0; 
otherwise  RATI01  is  set  equal  to  1. 
RATI01  is  calculated  by  dividing  the 
actual  transpiration  by  the  potential. 
RATI01  is  used  in  the  crop  growth  and 
root  growth  submodels  to  limit  dry  matter 
production  in  water  stress  situations. 

Figures  12.1,  12.2,  and  12.3  contain 
plots  of  predicted  transpiration, 
evaporation,  and  evapotranspiration  for  a 
typical  NTRM  model  run.   Note  the 
reduction  in  transpiration  from  about 
Julian  days  190  to  200.   This  was  caused 
by  a  temporary  drought.   The  crop  growth 
resumed  after  this  period  to  produce  the 
remaining  transpiration  pattern  for  the 
growing  season.   Figures  12.1,  12.2,  and 
12.3  illustrate  the  contributions  of 
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Figure   12.1. 

Transpiration  (corn)  from  an  NTRM  model 

run. 
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evaporation  and  transpiration  to  overall 
evapo transpiration.   Note  that  for  this 
particular  case,  the  general  shape  of  the 
ET  curve  was  determined  by  the 
transpiration  component.   The  shape  of 
the  evaporation  plot  was  determined 
mainly  by  the  precipitation  events  and 
the  surface  soil  hydraulic  properties. 
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Figure  12.3. 

Evapotranspiration  (corn)  from  an  NTRM 

model  run. 
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JULIAN  DAY  NUMBER 
Figure   12.2. 

Surface  evaporation   (corn)    from  an  NTRM 
model   run. 
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Chapter  13 
TILLAGE  SUBMODELS 

M.J.  Shaffer1 


Tillage  modifies  the  physical, 
biological,  and  chemical  properties  of 
the  soil,  which  in  turn  affect  crop 
yield.   The  manner  in  which  the  soil 
properties  are  modified  in  the  model 
takes  two  basic  forms.   First,  the  user 
can  specify  the  tillage  date,  depth,  and 
resulting  soil  properties  via  appropriate 
input  data,  or  secondly,  the  user 
specifies  only  the  tillage  date,  depth, 
and  type.   Tillage  submodels  are  then 
used  to  estimate  the  changes  in  existing 
soil  properties.   In  either  case,  once 
the  tillage  modifications  have  been  made, 
the  NTRM  model  uses  these  data  to  compute 
the  soil  water  balance  and  crop  yields. 

Changes  in  the  tilled  system  between 
tillage  events  are  also  considered. 
Again  users  have  two  options  available  to 
provide  this  information.   They  can  input 
them  as  baseline  conditions  for  the  soil 
properties  and  a  period  necessary  to 
reach  them  after  tillage,  or  they  can 
allow  the  tillage  submodel  to  compute 
these  changes . 

Table  13.1  contains  a  list  of  soil 
properties  that  the  model  directly 
changes  or  for  which  it  indirectly 
accepts  changes .   The  tillage  submodels 
use  regression  and  other  equations  to 
compute  changes  in  soil  bulk  density;  the 
soil  water  characteristic  curve; 
material  distributions;  and  surface 
properties  as  a  function  of  existing 
values  for  soil  texture,  soil  water 
content,  soil  organic  matter  content,  and 
surface  residues  for  various  tillage 
practices.   Computed  soil  bulk  density 
and  percent  organic  matter  are  then  used 
to  compute  the  updated  soil  water 
characteristic  curve,  soil  water  holding 


1  Research  soil  scientist, 
U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 


capacity,  and  soil  strength.   Carbon  and 
nitrogen  transformation  rates  are 
modified  as  a  function  of  tillage  depth. 

Mixing  of  soil  horizons  is  assumed  to  be 
either  complete  or  partial  depending  on 
the  tillage  practice  and  soil  conditions. 
The  same  holds  true  for  incorporation  of 
surface  residues.   The  resulting  horizons 
are  assumed  to  be  homogeneous  within 
their  own  boundaries.   This  means  that 
the  tillage  submodel  may  reconstruct  the 
location  of  horizon  boundaries  after  a 
tillage  event.   Soil  segment  boundaries 
used  for  internal  computations  are 
modified  if  the  soil  bulk  density 
changes . 

The  effect  of  tillage  on  soil  properties 
is  illustrated  in  figure  13.1.   In  this 
example,  certain  residues  are  partially 
incorporated  into  the  soil  and  two  soil 
horizons  are  completely  mixed.   The 
resulting  single  horizon  in  the  till  zone 
is  considered  to  be  homogeneous.   Some 
surface  residues  remain,  and  the  soil 
surface  properties  have  been  modified. 
The  deeper  soil  horizons  are  assumed  to 
be  unaffected,  although  deep  compaction 
could  also  be  considered. 


UNAFFECTED 


MODIFIED 
SURFACE 
PROPERTIES 


UNAFFECTED 


Figure  13.1. 

Modeling  tillage  effects  on  soil 

properties. 
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Table  13.1. 

Tillage-associated  soil  properties  modified  in  NTRM 

model 

Physical  properties    Biological  properties    Chemical  properties 


Soil  bulk  density 


Carbon  and  nitrogen 
transformation  rates 


Organic  matter 


Soil  water  charac- 
teristic curves 

Soil  water  holding 
capacity 


Microbial  mass  and 
distribution 

Root  mass  and 
distribution 


Residues 

N03"-N  and  NH4+-N 

Major  cations 
and  anions 


Soil  infiltration 
rate 


Cation  exchange  capacity 


Random  roughness 

Percentage  of 
residue  cover 


Solid  phase  salts 
Exchangeable  ions 


Soil  strength 

Soil  water 
content  (0) 
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Chapter  14 

SIMULATION  OF  INTERCEPTION,  SURFACE 

ROUGHNESS,  DEPRESSION  STORAGE,  AND  SOIL 

SETTLING 


D.R.  Linden  and  D.M.  Van  Doren,  Jr.1 


Introduction 


PlrM=°-08L-WC 


-CM 
where 


The  purpose  of  this  submodel  is  to 
hydrologically  distribute  rainfall  (or 
irrigation)  into  various  components  so 
that  the  net  addition  to  the  soil  water 
storage  system  can  be  obtained.   The 
model  is  developed  to  be  sensitive  to 
tillage  and  surface  residue  properties  of 
the  soil  system. 

During  rainfall,  precipitation  can  be 
partitioned  into  various  components  of 
water  balance.   The  following  section 
includes  a  discussion  of  this 
partitioning  into  the  interception  of 
rainfall  by  the  crop  canopy  and  residue 
mulches . 

Interception 

The  canopy  interception  of  rainfall 
subtracts  water  from  the  precipitation  so 
that  it  does  not  become  soil  water, 
drainage,  or  runoff.   This  water  is 
retained  by  the  crop  tissues  and  is  later 
evaporated.   Infiltration  and  canopy 
interception  can  occur  simultaneously 
since  a  growing  canopy  does  not 
completely  and  uniformly  cover  the  soil 
surface .   The  water  retained  by  the 
canopy  approaches  a  maximum  value  that 
represents  the  upper  limit  to  canopy 
interception.   The  upper  limit  to  canopy 
interception  given  by  Ritchie  (1972)  is 
expressed  as 


1  Linden  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108 
and  Van  Doren  is  a  professor 
emeritis,  Ohio  Agricultural 
Research  and  Development 
Center,  Ohio  State  University, 
Wooster,  OH  44691. 


W> 


XCM 


the  leaf  area  index, 

the  water  retained  in  the  canopy 

at  the  beginning  of  rainfall 

(cm) ,  and 

the  maximum  amount  of 

precipitation  that  is  retained  in 

the  canopy  (cm) . 


For  small  rains,  the  rainfall  amount  that 
is  retained  by  the  canopy  is  expressed  in 
the  same  proportion  as  the  amount  of  soil 
surface  coverage  as  manifested  in  a 
radiation  transmission  fraction  (Ritchie 
1972)  and  can  be  expressed  as 


PIc=(l-FT)P, 

f°r  P<[PICM/(1"FT)]' 
where 


[2 


■Vp  =  the  fraction  of  the  soil  surface 

area  not  shaded  by  the  canopy, 
P  =  the  rainfall  amount  (cm) ,  and 


ic 


the  canopy- intercepted  rainfall 
(cm)  . 


For  larger  rains 


P>PICM/(1"FT) 


the 


canopy  interception  would  be  limited  to 
the  maximum  water  that  the  canopy  can 
retain  (eq.  1) . 

Equation  1  relating  the  canopy 
interception  of  rainfall  to  the  leaf  area 
index  implies  that  all  types  of  crops 
behave  similarly  when  intercepting 
rainfall.   The  differences  in  canopy 
geometry  may  exist,  but  no  data  could  be 
found  to  distinguish  between  different 
types  of  crops. 

As  with  canopy  interception,  residue 
interception  of  rainfall  subtracts  water 
from  the  precipitation  so  that  it  does 
not  become  soil  water,  drainage,  or 
runoff.   This  water  is  retained  by  the 
crop-residue  tissues  and  is  later 
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evaporated.   Infiltration  and  mulch 
interception  can  occur  simultaneously 
since  the  residues  seldom  completely  and 
uniformly  cover  the  soil  surface.   The 
water  retained  by  crop-residue  tissue 
approaches  a  maximum  value  that 
represents  the  upper  limit  to  mulch 
interception.   The  upper  limit  to  mulch 
interception  given  by  Ritchie  (1972)  is 
expressed  as 


PlRM=0.08M-WRl 
where 


[3; 


M 


WV 


IRM 


the  dry  weight  of  residues  per 
unit  of  area  (metric  tons/ha) , 
the  water  retained  in  the  residues 
at  the  beginning  of  rainfall  (cm) , 
and 

the  maximum  amount  of 
precipitation  that  is  retained  in 
the  residue  mulch  (cm) . 


For  small  rains,  the  rainfall  amount  that 
is  retained  by  the  residues  is  expressed 
in  the  same  proportion  as  the  amount  of 
soil  surface  coverage  in  an  expression 
such  as 


PlR-(l-B)P, 


M 


f°r  P<[PlRM/(1-B)]' 


where 

B  =  the  fraction  of  the  soil  surface 
area  not  covered  by  residue 
material, 
P  =  the  rainfall  amount  (cm) ,  and 
Pj„  =  the  residue- intercepted  rainfall 
(cm)  . 

For  larger  rains  [p>Pj   /(l-B)j ,  the 

residue  interception  would  be  limited  to 
the  maximum  water  that  the  residues  can 
retain  (eq.  3) . 

Equation  3  relating  the  mulch 
interception  of  rainfall  to  the  dry 


weight  of  residues  implies  that  all  types 
of  crop  residues  behave  similarly  when 
intercepting  rainfall.   Differences 
between  stemmy  and  leafy  residues  may 
exist,  but  we  could  find  no  data  showing 
distinctions  between  different  types  of 
residues . 

Depression  Storage 

Depression  storage  and  associated  ponded 
infiltration  are  considered  after 
rainfall  ceases  and  some  water  has 
accumulated  on  the  surface  from  the 
infiltration  model.   In  this  model, 
depression  storage  is  the  maximum  amount 
of  water  that  can  be  retained  on  the 
surface  that  will  eventually  infiltrate. 
The  volume  of  space  on  the  surface 
available  for  storage  of  water  can  be 
computed  from  microrelief  data.   A 
reasonable  assumption  for  the  upper  plane 
of  this  storage  volume  is  the  equal  cut 
and  fill  or  the  equivalent  depth  plane. 
Further  consideration  of  the 
interconnections  of  ponded  depressions 
would  indicate  that  as  the  general  land 
slope  increases,  the  ability  of  the 
surface  to  retain  water  is  diminished.   A 
model  of  surface  depressions  (Linden 
1979)  was  used  in  developing  a  functional 
relationship  for  the  upper  limit  to 
depression  storage  as  a  function  of 
roughness,  RR,  and  the  general  land 
slope,  L.   This  model  analysis  did  not 
result  in  a  simple  functional 
relationship  but  can  be  expressed  in 
general  as 


D=f(RR,L) , 


5] 


where 

D 
RR 


the  upper  limit  to  depression 

storage , 

the  roughness  index  (standard 

deviation  of  height  measurements) 

and 

the  land  slope. 


91 


The  modeled  relationship  between  these 
three  variables  is  shown  in  figure  14.1. 
Depression  storage  has  a  maximum  value  of 
about  1  cm  and  decreases  as  roughness  and 
land  slope  increase. 

Point  surface  excess  water  (or  runoff)  is 
computed  as  the  surface  excess  generated 
from  the  infiltration  model  minus  the 
upper  limit  to  depression  storage. 

New  soil  properties  that  depend  on  the 
porosity  of  the  soil  are  computed  after 
calculation  of  updated  porosity.   It  is 
necessary  to  write  equations  for 
properties  as  functions  of  porosity  (or 
bulk  density)  to  accomplish  this  property 
updating. 

Soil  Surface  Smoothing 

Rough  surfaces  become  smoother  as  the 
soil  is  wetted.   Raindrop  impact  and 
microscale  erosion  and  soil  slumping 
account  for  this  smoothing.   This  model 
updates  the  roughness  condition  after 
each  rainfall  by  the  algorithm 
represented  by 

RR=RRI-(RRI-O.8)(l-EXP(-1.5B(cos0/A)E 

-0.015P))  [6] 

where 


RR  =  the  random  roughness, 

the  initial  random  roughness, 

the  exposed  surface  fraction, 

the  average  angle  of  inclination  of 

the  surface , 

the  surface  contact  area  per  unit 

of  horizontal  area, 

the  rainfall  energy,  and 

the  precipitation  amount. 


RRj  = 

B  = 

4  = 


A  = 


Equation  6  represents  an  exponential 
decrease  in  roughness  that  is  affected  by 
the  degree  of  protection  afforded  the 
surface  by  coverage  and  by  having 
inclined  soil  planes. 


Soil  Settling 

Freshly  tilled  soils  begin  settling  when 
exposed  to  wetting.   Structural  slumping 
is  probably  the  primary  mechanism  of  this 
settling.   An  algorithm  similar  to  the 
algorithm  for  smoothing  can  be  written  as 

a=aI-(aI-ac)(l-EXP(-1.5E-0.015P)) ,   [7] 

where 


a 

QI 

Qc 

E 

P 


the  porosity, 

the  initial  porosity, 

the  stable  porosity, 

the  rainfall  energy,  and 

the  rainfall  amount. 


It  should  be  noted  that  the  dependence  on 
rainfall  amount  and  rainfall  energy  is 
lower  than  the  previous  algorithm  (eq.  6) 
for  smoothing. 


Figure  14.1. 

The  upper  limit  to  depression  storage 
versus  general  land  slope  for  five  values 
of  random  roughness  (RR) . 
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Chapter  15 

CHEMICAL  EQUILIBRIA  SUBMODEL 

M.J.  Shaffer  and  S.C.  Gupta1 


Introduction 

Soil  system  reactions  that  are  more  rapid 
than  related  processes  can  be  simulated 
most  efficiently  using  a  chemical 
equilibria  approach.   In  this  situation, 
only  the  end  point  or  equilibrium  state 
need  be  considered,  eliminating  the  need 
to  follow  the  details  of  the  intermediate 
rate  processes.   This  technique  has  been 
applied  successfully  for  several  years  in 
soil  simulations  and  elsewhere.   The 
chemical  equilibria  submodel  currently 
used  in  the  NTRM  model  was  originally 
developed  by  Dutt  and  coworkers  in  the 
1960 's  (Dutt  1961,  Dutt  and  Tanj i  1962, 
Dutt  et  al .  1972)  and  has  seen  some 
refinement  and  extensive  application  in 
the  1970' s. 

This  submodel  uses  a  thermodynamic 
approach  to  predict  concentrations  of 
major  cations  and  anions  in  neutral- to - 
alkaline  soils.   Acid  soil  reactions  are 
not  considered.   Work  is  in  progress  to 
extend  this  submodel  into  the  acid 
region. 

Chemical  processes  being  simulated  in  the 
model  include  precipitation-dissolution 
of  CaC03  (lime)  and  CaS04«2H20  (gypsum); 
ion  pair  reactions  involving  CaS04° , 
MgS04° ,  NaS04~ ,  the  bicarbonate  buffer 


The  set  of  equilibrium  reactions  included 
in  the  model  can  be  written  as  follows: 


system;  and  ion  exchange  of  Ca 
Mg"1"4",  and  NH4+. 


Na^ 


1  Shaffer  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108  and 
Gupta  is  an  associate 
professor,  Soil  Science 
Department,  University  of 
Minnesota,  St.  Paul,  MN  55108. 


Gypsum  solubility- - 

Ksp=(Ca2+)(S042-)722, 
CaS04°  ion  pairing- - 

KD=((Ca2+)(S042-)/(CaS040))722. 
MgS04°  ion  pairing- - 

KD1=( (Mg2+) (S042' )/(MgS04° ) )722 

CaC03  solubility  and  bicarbonate  buffer 
system- - 


:i] 


3] 


K1PC02=(Ca2+)(HC03")272712, 
NaS04"  ion  pairing- - 

KD2=((Na+)(S042-)/(NaS04-))72,       I 
Ca-Na  ion  exchange 

J  aCa  /aNa=KCa-Na(ECa/ENa) - 
Ca-Mg  ion  exchange 

aCa/aMg=KCa-Mg(ECa/Emg) >  I 

Na-NH4  ion  exchange 

aNa/aNH4=KNa  -  NH4  <  %a/ENH4  )  ■  I 

Activity  coefficients,  7,  are  computed 
using  the  Debye-Huckel  equation  in  the 
form 

-log7i=0.509Zi2ru-/(l+.]~u~)  , 

where 

u=l/2XciZi2, 

Zj_   =  the  charge ,  and 
C^  =  the  concentration. 


[4] 


6] 


7] 


[9] 


10 
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The  solution  technique  used  to  solve  this 
set  of  nonlinear  simultaneous  equations 
is  as  follows: 

Each  equation  is  written  in  terms  of  its 
current  state  (not  at  equilibrium)  and 
stochiometric  changes  necessary  to 
achieve  equilibrium.   This  is  illustrated 
by  considering  equation  1  and  rewriting 
as , 


Ksp=((Ca2+)+AX)(S042_+AX)722 


[ii: 


This  equation  is  solved  for  AX  by  using 
the  Newton -Raphs on  technique  for 
polynomials , 


AXi+1-AXi-fi(X)/f(X) 


[12: 


The  completed  AX^+^  is  then  used  to 
update  the  appropriate  concentrations . 
The  other  equations  in  the  set  are  also 
solved  sequentially  using  the  same 
technique.   The  entire  procedure  is 
repeated  for  the  entire  equation  set 
until  each  AX  becomes  very  small  and  the 
system  approaches  equilibrium. 

In  practice,  the  Newton-Raphson  technique 
converges  very  rapidly,  and  a  minimal 
number  of  iterations  is  usually  needed  to 
achieve  satisfactory  convergence.   The 
chemistry  used  in  this  submodel  and 
additional  discussions  of  the  equation 
solution  techniques  are  reviewed  in 
Shaffer  and  Gupta  (1981). 


Model  Verification  and  Sensitivity 

This  submodel  has  had  extensive 
verification  and  use  over  the  past  two 
decades .   Numerous  practitioners  and 
research  people  from  the  United  States 
and  elsewhere  have  applied  some  version 
of  this  model  to  both  irrigated  and 
nonirrigated  agriculture.   The  U.S. 
Bureau  of  Reclamation  has  used  it  to 


study  environmental  impacts  of  projects 
involving  thousands  of  acres  and  millions 
of  dollars .   The  submodel  has  been 
adapted  to  simulate  the  equilibrium 
chemistry  of  sewage  effluents  applied  to 
soils  (Shaffer  and  Gupta  1981). 

Some  typical  comparisons  of  observed  and 
predicted  data  that  have  been  done  for 
this  submodel  are  illustrated  in  figures 
15.1  through  15.6.   Details  of  the 
studies  done  to  collect  these  data  and 
make  the  model  applications  can  be  found 
by  consulting  the  appropriate 
publications,  as  noted  in  the  figures. 
In  general,  the  submodel  has  a  reputation 
for  yielding  excellent  results  when  used 
within  its  limitations.   Several  model 
sensitivity  plots  are  shown  in  figures 
15.7-15.11  to  illustrate  the  response  of 
this  submodel  to  various  applied  inputs. 


Literature  Cited 

Dutt,  G.R.   1961.   Quality  of  percolating 
waters,  No.  1.   Development  of  a  computer 
program  for  calculating  the  ionic 
composition  of  percolating  waters. 
Contribution  No.  50,  Water  Resources 
Center,  Department  of  Irrigation, 
University  of  California,  Davis,  35  p. 

Dutt,  G.R.,  M.J.  Shaffer,  and  W.J.  Moore. 
1972.   Computer  simulation  model  of 
dynamic  bio-physiochemical  processes  in 
soils.   Univ.  Ariz.  Tucson  Agric.  Exp. 
Stn.  Tech.  Bull.  196,  101  p. 

Dutt,  G.R.  and  K.  Tanj i .   1962. 
Predicting  concentrations  of  solutes  in 
water  percolating  through  a  column  of 
soil.   J.  Geophys.  Res.   44:3437-3439. 

Shaffer,  M.J.  and  S.C.  Gupta.   1981. 
Hydrosalinity  models  and  field 
validation.   In  I.K.  Iskandar  (ed.), 
Modeling  Wastewater  Renovation,  Land 
Treatment,  p.  136-177,  John  Wiley  &  Sons, 
New  York. 


95 


No  (meq/Jt ) 
8  |  12 


Figure   15.1. 

Observed  (dashed  lines)  and  predicted 
(solid  lines)  sodium  distributions  after 
2  and  18  irrigations  with  water 
containing  11  meq/liter  sodium. 


Figure  15.2. 

Observed  (dashed  lines)  and  predicted 
(solid  lines)  chloride  distributions 
after  2  and  6  irrigations  with  water 
containing  18  meq/liter  chloride  during 
the  first  two  irrigations  and  21  meq/ 
liter  during  the  last  six. 
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Figure   15.3. 

Predicted  and  observed  selected  water 
quality  parameters,  Navajo,  WA. 
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Figure  15.4. 

Present  observed  and  predicted 
conditions,  South  Montezuma  Area,  CO. 
From  Shaffer  et  al .  (1977). 
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Figure   15.5. 

Observed  and  predicted  total  dissolved 
solids  (TDS)  and  major  constituent 
concentrations  in  drainage  effluent, 
Suisun  March  study  site,  CA. 
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Figure  15.6. 

Calculated  vs .  measured  exchange  percent 

for  Ca"1^",  Mg++,  and  Na+. 
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Figure  15.7. 

Total  dissolved  solids  (TDS) 

concentration  vs.  depth  over  time. 


97 


0.0 


■~- 

50.0 

c 

U 

*-■ 

r 

i- 

0. 

H 

10O.U 

iso.o 


«0.g 

^  ["6/h  3QD.Q 


Figure   15.8. 

Sulfate  concentrations  vs.  depth  over 

time. 
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Figure  15.10. 

(a)  Breakthrough  curves  for  sodium 
adsorption  ratio  (SAR)  ,  and  (b)  changes 
in  exchangeable  sodium  percentage 
profiles  influenced  by  soil  type. 
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Figure   15.9. 

Changes  in  exchangeable  sodium  percentage 
profiles  and  breakthrough  curve  for 
sodium  adsorption  ratio  (SAR)--loam  soil 
irrigated  with  various  waters. 
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Figure  15.11. 

Changes  in  sodium  adsorption  ratio  (SAR) 

(a)  and  exchangeable  sodium  percentage 

(b)  profiles  as  influenced  by  depth  of 
reactive  profile. 
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Chapter  16 

SOLUTE  TRANSPORT  SUBMODEL 

M.J.  Shaffer1 


Introduction  and  Theory 

Transport  of  soluble  species  in  the 
unsaturated  zone  is  handled  in  the  NTRM 
model  via  the  submodel  FL1 .   This  route 
uses  data  on  soil  water  flux  and  soil 
water  content  provided  from  the 
unsaturated  flow  submodel  or  elsewhere  to 
calculate  solute  fluxes  in  and  out  of 
each  soil  segment  for  each  time  step. 
The  algorithms  used  can  be  summarized  as 


J   t,J-l  JJ-1     t  „J  . 

AmT .    —C.  • F    - C . • F  ,  and 

11  l 

[for  positive  flux] 

AmT^=Ct'J.FJ-Ct'J+1.FJ+1 
li       l 


[i: 


[2 


for  negative  flux] 


where 

AmT  = 

C  = 

F 
t 
J 


the  change  in  solute  mass  for 

segment  J , 

concentration  of  solute 

species  "i", 

water  flux, 

the  time  step  superscript,  and 

the  soil  segment  superscript. 


Equations  1  and  2  are  written  for  general 
flow  in  the  downward  or  upward 
directions.   The  FL  submodel  also  allows 
simultaneous  upward  and  downward  flow 
into  any  soil  segment.   This  gives  more 
flexibility  in  the  solute  redistribution 
process . 


so-called  "numerical  dispersion"  into  the 
simulation.   This  allows  a  first 
approximation  of  dispersion  without  the 
circular  procedure  of  removing  the 
numerical  effects  and  then  re -adding 
dispersion  via  a  Ficks  Law  dispersion 
term. 

Studies  have  shown  that  the  first 
approximation  is  satisfactory  for  most 
soils.   The  amount  of  dispersion  being 
introduced  can  be  controlled  to  some 
extent  by  the  number  of  soil  segments 
used  in  the  profile.   If  necessary,  it 
would  be  possible  to  include  a  circular 
dispersion  algorithm  in  the  NTRM  model. 


Submodel  Verification  and  Sensitivity 

Figures  16.1,  16.2,  and  16.3  contain 
predicted  and  observed  data  for  Na+,  CI", 
and  NO3",  respectively,  using  the  solute 
transport  algorithm  in  the  NTRM  model. 
The  agreement  with  experiment  shown  in 
these  graphics  is  typical  of  results  that 
can  be  expected  from  the  model.   Note 
that  the  curves  shown  in  figures  16 . 1  and 
16.2  represent  results  obtained  after 
both  short-term  and  longer  term  leaching. 
The  N03_-N  data  displayed  in  figure  16.3 
were  collected  at  the  end  of  a  growing 
season,  and  the  model  simulation  run  was 
for  the  entire  season. 


Solute  transport  due  to  dispersion  or 
diffusion  is  not  handled  directly  in  the 
model.   Instead,  the  use  of  finite  soil 
segments  with  complete  mixing  assumed  at 
each  step  in  time  and  space  introduces 

1  Research  soil  scientist, 
U.S.  Department  of 
Agriculture ,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108. 
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Figure  15.1. 

Observed  (dashed  lines)  and  predicted 
(solid  lines)  sodium  distributions  after 
2  and  18  irrigations  with  water 
containing  11  meq/liter  sodium. 
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Figure  15.2. 

Observed  (dashed  lines)  and  predicted 
(solid  lines)  chloride  distributions 
after  2  and  6  irrigations  with  water 
containing  18  meq/liter  chloride  during 
the  first  two  irrigations  and  21  meq/ 
liter  during  the  last  six. 


Figure  16.1. 

Observed  (dashed  lines)  and  predicted 
(solid  lines)  sodium  distributions  after 
2  and  18  irrigations  with  water 
containing  11  meq/liter  sodium. 


101 


Chapter  17 

CROP  RESIDUE  SUBMODEL 

M.J.  Shaffer  and  J.A.E.  Molina1 


Introduction 

Simulation  of  crop-residue  interactions 
with  soil  properties  and  crop  yield  occur 
in  several  places  in  the  NTRM  model. 
Crop -residue  placement  and  incorporation 
are  illustrated  in  figure  17.1.   Note 
that  residues  are  allowed  on  the  soil 
surface  as  well  as  being  incorporated 
into  the  soil  material.   Processes 
modeled  with  the  surface  residues  include 
(1)  nitrogen  and  carbon  transformations 
as  functions  of  residue  temperature  and 
water  content  and  (2)  leaching  of  soluble 
residue  constituents. 

Concentration  units  used  with  the  surface 
residues  are  ng   residue  material/cm^-  of 
soil  surface.   This  is  in  contrast  to  the 
subsurface  units  of  /xg  residue  material/ 
Hg   soil. 

When  the  residues  are  incorporated  into 
the  soil,  the  model  considers  nitrogen 
and  carbon  transformations;  leaching  and 
redistribution  of  soluble  constituents; 
and  the  effects  of  residues  on  physical, 
biological,  and  chemical  properties  of 
the  soil. 

Surface  residue  water  contents  and 
temperatures  are  computed  in  the  surface 
residue  submodel.   Residue  carbon  and 
nitrogen  transformations  are  calculated 
in  the  nitrogen  transformation  submodel. 
Leaching  and  redistribution  of  soluble 
residue  constituents  are  computed  by  the 
solute  transport  submodel.   The  tillage 
submodels  handle  the  effects  of  residue 


1  Shaffer  is  a  research  soil 
scientist,  U.S.  Department  of 
Agriculture,  Agricultural 
Research  Service,  Soil  and 
Water  Management  Research 
Unit,  St.  Paul,  MN  55108  and 
Molina  is  a  professor,  Soil 
Science  Department,  University 
of  Minnesota,  St.  Paul,  MN 
55108. 


incorporation  on  soil  properties.   Root 
slough  produced  in  the  root  growth 
submodel  becomes  residue  material 
available  for  decay  and  for  generation  of 
miscellaneous  effects  on  soil  properties. 


Surface  Residues 

Carbon  and  nitrogen  transformation 
(decay)  processes  involving  residues  on 
the  soil  surface  are  modeled  using  the 
nitrogen  transformation  submodels 
(chapter  7,  "Carbon  and  Nitrogen 
Transformations  in  Soil,  Submodel 
NCSOIL"),  and  the  transformation 
submodel,  TRNSFM,  reported  by  Shaffer  et 
al.  (1977).   Rate  coefficients  for  the 
transformation  equations  are  specific  for 
residues  located  on  the  surface.   Note 
that  the  units  and  values  of  these 
surface  coefficients  are  slightly 
different  from  their  subsurface 
counterparts .   These  rates  are  taken  as  a 
function  of  residue  water  content  and 
residue  temperature.   Air  temperature, 
soil  temperature,  and  residue  type  are 
used  to  compute  a  representative 
temperature  for  the  residues  using  the 
relationship 

Tr=ABS(Tair-Tsoil)/2.Kres.Kth.Kcov     [1] 

where 


T 

r 

—  temperature  of  residue 

(°C), 

T  • 
-"-air 

=  air  temperature  just  above 
residue  (°C) , 

^soil 
res 

=  soil  surface  temperature 

(°C), 
=  residue  coefficient, 

Kth 

=  residue  thickness 

coefficient,  and 

v 

^COV 

=  residue  coverage 
coefficient. 
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Figure   17.1. 

Processes  modeled   involving  crop 

residues. 


residue  temperature  and  water  content  are 
used  by  the  nitrogen  transformation 
submodels.   NH^-N  and  N03~-N  mineralized 
during  residue  decay  are  available  for 
both  leaching  into  the  soil  profile  and 
removal  with  surface  runoff. 


Incorporated  Residues 

Nitrogen  and  carbon  transformations 
involving  residues  incorporated  into  the 
soil  are  computed  using  the  nitrogen 
transformation  submodels,  (chapter  7, 
"Carbon  and  Nitrogen  Transformations  in 
Soil,  Submodel  NCSOIL"),  submodel  TRNSFM, 
and  specific  rate  coefficients.   Soil 
temperatures  and  soil  water  contents  used 
in  these  calculations  are  obtained  from 
the  soil  temperature  and  unsaturated  flow 
submodels . 


Residue  water  content  is  based  on  a  water 
budget  for  the  residue,  using  a  residue 
field  capacity  after  a  precipitation 
event  and  the  evaporation  rate  for  the 
residue  surface.   These  values  for 


The  effects  of  incorporated  residues  on 
soil  properties  such  as  bulk  density, 
surface  roughness,  and  surface 
reflectancy  are  computed  using  regression 
and  other  equations . 
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